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Abstract 


During  the  grant  period,  we  have  developed  the  high  order  methods,  namely  spectral 
methods,  WENO-Z  finite  difference  methods  and  the  Hybrid  spectral- WENO  finite 
difference  scheme  for  supersonic  reactive  and  non-reactive  flows.  High  order  methods 
are  needed  in  these  kind  of  problems  because  the  numerical  solutions  contain  small 
scales  that  are  important  in  the  simulations  of  turbulence  in  interface  instability  and 
combustion  in  modeling  combustors.  Also  since  most  of  the  relevant  problems  are  not 
steady,  long  term  integrations  have  to  be  carried  out.  It  is  well  known  that  high  order 
accuracy  methods  are  needed  in  this  situation. 


We  consider  the  problem  of  supersonic  reactive  flow  in  recessed  cavities  and  simulate 
it  by  spectral  multidomain  techniques.  This  is  important  as  a  flameholder  mechanism 
in  scramjets.  We  have  developed  a  multi-domain  spectral  method  with  stable  and 
conservative  penalty  interface  conditions  for  the  numerical  simulation  of  supersonic 
reactive  recessed  cavity  flows  with  homogeneous  grid. 


The  multi-domain  hybrid  Spectral- WENO  finite  difference  method  is  introduced  for 
the  numerical  solution  of  two  dimensional  nonlinear  hyperbolic  systems  in  a  Cartesian 
physical  domain  which  is  partitioned  into  a  grid  of  rectangular  subdomains.  The 
main  idea  of  the  Hybrid  scheme  is  to  conjugate  the  spectral  and  WENO  methods 
for  solving  problems  with  shock  or  high  gradients  such  that  the  scheme  adapts  its 
solver  spatially  and  temporally  depending  on  the  smoothness  of  the  solution  in  a 
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given  subdomain.  Built  as  a  multi-domain  method,  an  adaptive  algorithm  is  used  to 
keep  the  solutions  parts  exhibiting  high  gradients  and  discontinuities  always  inside 
WENO  subdomains  while  the  smooth  parts  of  the  solution  are  kept  inside  a  spectral 
one,  avoiding  oscillations  related  to  the  well-known  Gibbs  phenomenon  and  increasing 
the  numerical  efficiency  of  the  overall  scheme.  A  higher  order  version  of  the  multi¬ 
resolution  analysis  proposed  by  Harten  is  used  to  determine  the  smoothness  of  the 
solution  in  each  subdomain.  The  Hybrid  method  is  applied  to  the  two-dimensional 
Shock- Vortex  Interaction  and  the  Richtmyer- Meshkov  Instability  (RMI)  problems. 


We  also  have  developed  a  different  approach  in  attempting  to  get  more  meaningful 
results  is  to  model  statistically  those  scales  that  can  not  be  resolved.  Methods  for 
modeling  those  scales  are  being  developed  and  applied.  Numerical  issues,  as  boundary 
conditions,  are  being  addressed  and  new  boundary  procedure  is  being  presented. 


1  A  weighted  multi-domain  spectral  penalty  method  with  inhomogeneous 
grid  for  supersonic  injective  cavity  flows 


Spectral  methods  have  been  actively  used  in  the  computational  fluid  dynamics  com¬ 
munity  in  the  last  decades  due  to  the  merit  of  high  order  accuracy  maintained  for 
long  time  integration.  Spectral  methods  also  have  been  applied  to  highly  complex 
fluid  systems  and  have  been  proved  to  yield  accurate  solutions  even  with  the  stiff  or 
discontinuous  spatial  gradients.  These  systems  include  the  supersonic  shock  bubble 
interactions  [12],  the  supersonic  cavity  flows  [11],  and  etc.  The  difficulty  of  implement¬ 
ing  the  spectral  method  to  these  complex  fluid  systems  is  to  deal  with  the  stiff  or 
discontinuous  spatial  gradients  successfully.  The  discontinuous  solution  is  commonly 
found  in  most  high  speed  fluid  mechanical  systems.  The  spectral  approximation  of 
such  solutions  yields  spurious  oscillations  near  the  discontinuity,  known  as  the  Gibbs 
phenomenon.  These  Gibbs  oscillations  destroy  both  the  accuracy  and  stability  in  gen¬ 
eral.  The  essential  methodology  to  deal  with  such  oscillations  in  the  spectral  solution 
is  the  spectral  viscosity  or  Altering  methods  [5,16,20,23].  The  Altering  which  is  math¬ 
ematically  equivalent  to  the  spectral  viscosity  method  but  practically  more  efficient, 
is  used  to  stabilize  the  Aow  Aelds  over  the  time  integration.  The  Altering  reduces  the 
high  order  oscillations  by  attenuating  the  high  modes  in  the  solution  in  a  smooth 
manner.  The  Altering  method  can  be  applied  either  globally  or  locally.  By  applying 
the  Altering  locally  more  accurate  solution  is  obtained  in  the  smooth  region  and  the 
spurious  oscillations  in  the  neighborhood  of  the  discontinuity  is  considerably  reduced. 
Thus  it  is  desirable  to  separate  the  locally  non-smooth  regions  from  the  global  smooth 
region.  A  multi-domain  spectral  method  has  been  developed  to  address  this  problem 
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[8,11,13,14,17-19].  Once  the  physical  domain  is  split  into  multiple  subdomains,  the 
proper  interface  conditions  should  be  imposed  at  the  domain  interfaces.  The  simplest 
condition  is  obtained  using  the  averaging  method.  With  the  averaging  method,  the 
flow  field  at  the  domain  interface  is  given  by  the  average  of  the  two  adjacent  solu¬ 
tions  across  the  interface.  Thus  the  continuity  of  the  solution  is  ensured  with  the 
averaging  method.  Although  this  method  is  simple  and  efficient  to  be  implemented, 
it  may  cause  the  generation  of  nonphysical  solutions  at  the  interface  if  the  two  ad¬ 
jacent  subdomains  have  different  order  of  approximation  or  different  domain  length, 
i.e.  if  the  grid  system  is  inhomogeneous.  We  define  the  grid  inhomogeneity  as  the  grid 
configuration  such  that  the  grid  resolutions  between  the  adjacent  subdomains  across 
the  domain  interface  are  different.  Such  difference  can  be  obtained  by  having  each 
domain  have  different  order  of  polynomials  or  different  physical  domain  length.  If  the 
distribution  of  the  subdomain  grids  is  inhomogeneous,  the  stable  interface  conditions 
derived  for  the  homogeneous  grid  system  are  not  enough  and  one  needs  to  find  the 
conditions  with  which  the  spatial  inhomogeneity  can  be  addressed  properly. 


At  the  domain  interface  of  two  adjacent  subdomains  which  have  the  degree  of  poly¬ 
nomials,  N\  and  N2  in  x— direction,  the  ratio  of  the  grid  spacing  Ax±  to  Ax2  is 
approximately  given  by 

Ax2  =  An  Nf 

A.r,  A7  '  Xf  1  j 

where  A1  and  A11  are  the  domain  lengths  of  the  two  subdomains.  If  the  grid  spacing 
ratio  ^  is  different  and  far  from  unity,  we  consider  it  as  the  inhomogeneous  grid.  If 
Ai^.  =  1,  i.e.  if  the  grid  is  homogeneous,  the  averaging  method  can  play  an  efficient 
role  as  a  stable  interface  condition.  However,  if  the  ratio  is  far  from  unity,  the  simple 
averaging  interface  condition  can  cause  the  growth  at  the  domain  interface.  In  the  real 
computation,  the  values  of  N\ ,  N2  and  A7,  A77  are  chosen  such  that  becomes  close 
to  unity  but  in  general  ^  1.  The  current  work  centers  around  the  development 
of  the  method  dealing  with  the  solution  in  inhomogeneous  grid,  i.e.  when  ^  1. 


In  [11],  we  performed  a  2D  direct  numerical  simulation  (DNS)  of  the  recessed  cavity 
with  the  multi-domain  spectral  penalty  method  under  the  condition  that  the  grid  is 
homogeneous.  For  the  homogeneous  grid  system,  the  size  and  the  number  of  collo¬ 
cation  points  in  each  subdomain  are  the  same  in  each  dimension.  In  this  study,  we 
extend  the  previous  work  to  the  inhomogeneous  grid  system  to  consider  the  injector- 
cavity  system  with  the  local  hydrogen  fuel  injector.  The  crucial  part  of  the  DNS  of 
the  injector-cavity  system  is  to  resolve  the  hydrogen  jet  injector  without  causing  any 
instability  or  nonphysical  growing  modes  at  the  domain  interfaces.  The  ratio  of  the 
injector  to  the  cavity  length  scale  is  about  (^(lCr1).  We  use  a  smaller  subdomain 
with  higher  order  polynomials  to  resolve  the  narrow  jet.  In  Figure  1  the  local  do- 


3 


main  configuration  is  given  for  the  cavity  flameholder  with  (left  figure)  and  without 
(right  figure)  the  injector.  The  local  domain  configuration  shown  in  the  right  figure 
is  the  typical  domain  system  used  in  [11]  for  which  grid  system  is  homogeneous.  The 
grid  system  in  the  left  figure  is  inhomogeneous  as  the  local  injector  is  in  the  narrow 
domain.  In  [11],  the  stability  analysis  has  been  done  with  the  assumption  that  each 
subdomain  has  the  same  length  but  can  have  different  polynomial  orders.  With  differ¬ 
ent  polynomial  orders,  the  stability  is  still  maintained.  In  this  work,  we  further  show 
that  the  stability  can  be  also  maintained  with  the  different  domain  length.  For  the 
2D  DNS,  we  use  the  grid  inhomogeneity  mainly  due  to  the  different  domain  length 
with  the  same  polynomial  order  in  each  domain. 


The  stability  conditions  can  be  derived  with  the  inhomogeneous  grid,  but  these  are 
only  necessary  conditions.  Consequently  there  can  be  a  non-reflecting  mode  at  the  do¬ 
main  interface  which  can  yield  a  growth  in  time.  A  weighted  spectral  penalty  method 
is  proposed  in  order  to  minimize  such  nonphysical  reflecting  modes  at  the  inhomo¬ 
geneous  grid  interfaces.  We  note  that  in  [8]  the  multidomain  spectral  method  has 
been  also  used  for  the  localized  incompressible  stratified  turbulence  flows  in  which  the 
strong  adaptive  averaging  method  has  been  used  with  the  spectral  filtering  technique. 
In  our  work,  we  observe  that  the  averaging  method  with  the  filtering  technique  does 
not  yield  a  smooth  profile  across  the  inhomogeneous  grid  interface  for  the  compress¬ 
ible  supersonic  flows.  This  is  due  to  the  fact  that  the  averaging  does  not  guarantee  the 
stability  for  the  inhomogeneous  grid  as  shown  in  the  following  section.  The  character¬ 
istic  decomposition  type  interface  conditions  also  fail  to  provide  the  smooth  solution 
across  the  inhomogeneous  grid  interfaces.  The  major  developments  of  the  current 
work  are  as  follows.  1)  A  generalized  conservative  and  stable  penalty  conditions  are 
derived  for  the  inhomogeneous  grid.  2)  The  weighted  spectral  penalty  method  is  de¬ 
veloped  to  minimize  the  non-physical  growth  modes  at  the  inhomogeneous  domain 
interfaces. 


Using  the  proposed  weighted  penalty  method,  we  carry  out  the  2D  DNS  of  the 
injector-cavity  system  with  various  recessed  angles.  Cavity  is  an  efficient  flame-holder 
of  scramjet  engine  as  it  generates  the  self-sustained  recirculation  region.  The  hot 
radicals  from  the  chemical  reactions  residing  in  the  recirculation  region  reduce  the 
induction  time  and  consequently  maintain  the  auto-ignition.  For  the  continuous  auto- 
ignition  and  better  fuel  efficiency,  such  recirculation  region  should  be  stable  for  long 
time.  In  addition  to  the  recirculation,  the  self-sustained  acoustic  oscillations  bouncing 
back  and  forth  inside  cavity  disturb  the  recirculation  generating  pressure  fluctuations. 
The  geometry  of  cavity  is  an  important  parameter  for  maintaining  a  stable  recircu¬ 
lation  while  reducing  the  pressure  oscillations.  It  is  shown  in  [11]  that  the  recessed 
cavity  flameholder  reduces  the  pressure  fluctuations  inside  cavity  more  considerably 
than  the  normal  wall  cavity.  In  this  research  we  verify  qualitatively  that  the  recessed 
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cavity  increases  the  stability  of  the  recirculation  and  reduces  the  pressure  fluctuations 
inside  the  recessed  cavity  with  the  hydrogen  injector. 


Air 

M=1 .91 

Re=3.95E7  (1/m) 

Pt=2.828522(atm) 

Tt=830.6(K) 

L/D=4 

::::> 
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Fig.  1.  Left:  Local  domain  configuration  of  the  normal  injector-cavity  flame-holder.  Right: 
Local  domain  configuration  of  the  normal  cavity  without  injector.  The  initial  physical  con¬ 
figuration  of  the  injector-cavity  flame-holder  is  given  in  the  legend  box,  where  M  denotes 
the  Mach  number,  Re  the  normalized  Reynolds  number,  Pt  the  baseline  total  pressure,  Tj 
the  baseline  total  temperature,  and  L/D  the  length  to  depth  ratio  of  cavity. 

The  paper  is  organized  as  follows.  In  Section  2,  a  Legendre  multi-domain  spectral 
method  with  inhomogeneous  grid  is  explained.  Stability  and  conservativity  are  de¬ 
rived  with  the  grid  inhomogeneity.  The  generalized  penalty  interface  conditions  are 
derived  accordingly.  The  weighted  penalty  method  is  proposed  and  various  examples 
are  illustrated.  In  Section  3,  the  governing  equations  and  the  injector-cavity  system 
are  briefly  described.  The  numerical  results  from  the  simulation  of  the  supersonic  cav¬ 
ity  flame-holder  are  provided.  In  this  section,  we  verify  that  the  proposed  weighted 
spectral  penalty  method  is  applied  successfully  with  the  inhomogeneous  grid  of  the 
injector-cavity  system.  The  pressure  fluctuations  both  in  the  normal  and  30°  wall 
cavities  are  presented.  Concluding  remarks  with  a  brief  future  work  outline  are  given 
in  Section  4. 


2  Multi-domain  spectral  method  with  inhomogeneous  local  mesh  refine¬ 
ment 


2.1  Conservative  spectral  penalty  7nethods  for  inhoinogeneous  grid 


In  this  section,  we  will  consider  the  spectral  penalty  method  for  the  inhomogeneous 
grid.  We  first  consider  the  following  one-dimensional  conservation  laws: 


dq(x,t )  df(q(x,t )) 

v  ;  +  v '  =  0,  x  e  R,  t  >  0. 


dt 


dx 


(2) 


Here  q(x,t )  is  the  state  vector  and  f(q(x,t ))  is  the  flux  vector.  In  [6]  the  conservative 
multi-domain  Legendre  method  was  proposed  to  approximate  (2)  on  the  Legendre 
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Gauss- Lobatto  collocation  points  and  such  formulation  was  successfully  used  in  [11], 
In  [6],  theorems  for  the  multi-domain  Legendre  penalty  method  have  been  provided 
under  the  assumption  that  each  subdomain  has  the  same  domain  length  but  the 
polynomial  orders  of  approximation  can  be  different.  The  different  polynomial  order  in 
each  domain  makes  the  grid  system  inhomogeneous.  Thus  the  same  polynomial  order 
has  been  used  in  [11]  to  avoid  any  numerical  artifacts  due  to  such  inhomogeneity.  In 
this  paper,  we  will  further  generalize  the  previous  formulation  for  fully  inhomogeneous 
grid  system. 


For  simplicity,  we  consider  two  subdomains  tt1  =  [xx,,0]  and  Q77  =  [0,  xr],  for  which 
the  domain  interface  is  at  x  =  0.  In  [11],  xl  =  —xr  =  —2,  but  xl  and  xr  can  be 
different  in  this  paper.  Furthermore  the  left  domain  uses  the  polynomial  order  of  N 
and  the  right  domain  of  M  and  N  is  not  necessarily  the  same  as  M.  The  Legendre 
multi-domain  spectral  penalty  method  is  then  given  by 


do1  dl1  fia1) 

m  +  d,  =  *))  +  svwN)  + 

n<3^(V)[/+te( 0,t))  -  ./+(g£(0,t))]  + 

r2QN(x)[f-(qIN(0,t))  -  0,t))], 

da11  dl11  fia11) 

-Jf  +  Mq^m)  =  B(ti(xR,  t))  +  SV(q“)  + 

T3QM(x)[f+{q^(0,t))  -  /+(g^(0,t))]  + 
uQM(x)[r{q^(0,t))  -  f~(qIN{0,  £))].  (3) 


Here  q!N  denotes  the  numerical  approximation  of  q{x,  t )  in  Legendre  polynomial  of 
order  N  in  IT  and  q1^  of  order  M  in  Qn .  B  is  the  boundary  operator  at  the  end 
points,  i.e.  x  =  xl,  xr  and  SV  is  the  spectral  vanishing-viscosity  terms.  I'N  and  I'J  are 
the  Legendre  interpolation  operators  for  the  left  and  right  subdomains  respectively. 
Q N  and  QM  are  the  polynomials  of  order  N  and  M  respectively  defined  to  vanish 
at  the  collocation  points  except  at  the  boundary  or  interface  points,  that  is,  for  Q1 , 
Qn^xi)  —  0  for  i  —  1,  •  •  •  ,  N  —  1  and  Qn(xi)  =  1  for  i  —  0,N.  The  positive  and 
negative  fluxes  f+  and  f~  are  defined  by 

/±  =  /  SA±S~1dq,  (4) 


with 

A=^-  =  SAS-1.  (5) 

dq 

The  Jacobian  matrix  A  is  assumed  to  be  symmetric.  A+  and  A~  are  the  diagonal 
matrices  composed  of  positive  and  negative  eigenvalues  of  A  respectively  such  as 


6 


A  =  A+  +  A~.  S  and  A  are  the  variables  related  to  the  characteristics  and  its  direction 
of  propagation.  t1,T2,t3  and  r4  are  the  penalty  parameters  and  all  are  constants.  As 
in  [11],  we  assume  that  the  boundary  terms  and  the  spectral  vanishing- viscosity  terms 
do  not  cause  any  instabilities  and  they  do  not  appear  in  the  following  analysis.  For  the 
following  theorems  we  define  the  discrete  Legendre  norm  (p,  q )  N  :=  YliLo 
Xi  are  the  Legendre  Gauss- Lobatto  collocation  points  and  cu*  =  v(Ar+i)|Xjv(g(3;-))]2  w^iere 
£  is  the  map  from  x  to  the  Legendre  Gauss-Lobatto  points  over  [—1, 1],  If  pq  G  P2n- i, 
the  discrete  sum  is  exact,  i.e.  (p,  q)^  =  f\  p(£(x))q(£(x))d£.  In  the  following  analysis, 
we  define  the  weight  vector  u VN  as  the  weight  vector  in  with  IV  +  1  components 
such  as  =  (oUq,  ■  ■  ■  ,uIN)T.  We  note  that  without  the  vector  symbol  denotes 
the  last  component  of 


The  scheme  given  in  (3)  is  conservative  if  xl  =  —  Xr  =  —2  and  the  penalty  parameters 
satisfy  the  following  conditions 


-  T3U%  =  1,  72 Cod  -  r4u =  1.  (6) 

Here  we  note  that  there  is  a  typo  in  the  first  equation  of  Eq.  (23)  on  page  332  in 
[11].  The  second  term  in  the  equation  should  not  be  r4u;{|  but  should  be  T3cy[(  as 
given  in  the  theorem  above. 


The  scheme  (3)  is  stable  if  xl  —  —xr  =  —2  and  the  penalty  parameters  satisfy  the 
followings 


2t4c4  <  1,  2t2uin  >  1,  2t3cu"  <  -1,  2 t4u;"  >  -1,  (7) 

nujrN  -  =  1,  r2o;^  -  r4u;"  =  1.  (8) 

Here  note  that  the  scheme  is  stable  even  though  the  grid  system  is  inhomogeneous, 
i.e.  N  ±  M. 


If  each  subdomain  has  the  same  domain  interval  A,  then  the  stability  conditions  are 
given  by,  defining  A2  =  2/ A, 


2r4u;^  <  A2,  2r2cod  >  A2,  2 t3lo“  <  -A2,  2r4o;"  >  -A2, 

Ti ujf  -  r3u^  =  A2,  t2(xin  -  r4u;"  =  A2. 
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Proof:  The  proof  is  done  easily  using  the  fact  that  for  Q1 , 


since  J^/(g^r)  G  P2N-1  and  in  the  same  way  for  QH . 


f(/(o)-/(mo9) 


If  the  interval  of  each  subdomain  is  different,  then  the  scheme  (3)  is  conservative  if 
the  following  conditions  are  satisfied. 


^  ,  A"  „ 

yh uN  -  — r3 ujm 


1, 


A1 

- ToUJ 

2 


1 

N 


A11 

2 


,  JI 

A  ^  M 


1. 


(10) 


Proof:  Multiply  the  equations  for  qrN  and  q"  in  (3)  by  n VM  =  (uq,  ■  ■  ■  )x>IN)T  and 
c d1^  =  (cUq1  ,  •  •  •  ,  Um)T-  Then  using  the  Legendre  quadrature  rule  we  have 


r°  JK 

JxL  dt 


dx  + 


fXR  Qq 


dt 


dx 


+ 

+ 


df"rJr  -  H  dJ™- 

JxL  3x  '  Jo  dx 

y-n[/+te(0,t)V;  - 

y-rs[/+(gM(o,t)V"- 

Ar4[/V9"(0,«)K'- 


/+(tf(0,t)K] 
/“(?«( M))<4] 


Using  the  fact  that  /  ^-dx  =  f+  +  f  and  uJq  =  cn",  we  have  the  RHS  of  the  above 
equation  without  the  boundary  terms  become 


A7  A11  A11 

RHS  =  f+(qIN(0,t))[—T1uIN  -  —t3u%  -  1]  +  /+(9£(0,t))[— t3u“ 

/”te(0,t))[y-A<  -  yj-T4o;"  -  1]  +  /'(^"(O,  t))[^-r4^ 

For  any  /±(0,t)  (note  that  f±(qIN(0,t ))  7^  f±(q^(0,t))  in  general),  the  RHS  vanishes 
if  the  conditions  Eq.  (10)  are  satisfied.  □ 


A'  , 

— yu^jv  + 

A7  , 
- —t2^n  + 


The  scheme  (3)  is  stable  if 


2  2  2  2 
2ti<4  <  — ,  2t30 <  —  ^77’  2t2(4  —  2t40 4  >  —  ^77,  (H) 

2  2  2  2 

-  r3u$)2  -  2(tiu4^77  -  t3u$— 7)  +  <  0, 

(t2^n  ~  h4)2  -  2(^4^  -  +  ^7^77  ^  °>  (12) 


Proof:  Multiplying  each  of  Eq.  (3)  by  g7  and  g77>  then  the  energy  E(t)  =  -p 7  g2(x,  t)c7r+ 

-pr  /qR  q2(x,  t)dx  becomes 


ld£(t)  ,  j  2  1.  ,  .  j  n.  ,  .  n  2  L  , 

2  ^  =  Mi^iv  -  ^7 2^ao  -  (ri^7v  +  t's^mJTo  +  (a^m  +  ^j7-)/?o 

2  1  2  1 
+  (T2<4  -  ^7  2)^0  -  (^<4  +  rAuj^)%  +  (r44r  +  ~£jj~j)Po  , 

where  a*  =  ((4(0,  t))T7l±g^(0,  t),  $  =  (4(0,  t))TA±q^(0,  t),  and  7^  =  (4(0,  t))r7l±g^(0,  t). 
To  make  the  RHS  less  than  or  equal  to  zero,  first,  the  coefficients  of  the  2nd  order 
terms  corresponding  the  positive  flux  (negative  flux)  should  be  non-positive  (non¬ 
negative)  which  provides  the  first  conditions  in  Eq.  (12).  Also,  the  determinant  each 
of  the  quadratic  equations  should  be  non-positive.  This  provides  the  second  conditions 
Eq.  (12).  □ 


Fig.  2.  Stability  regions.  Left:  A7  =  A77  =  2.  Right:  A7  =  2  and  A77  — >  00. 

Figure  2  shows  the  stability  regions  for  TiU^  and  r2 with  various  A77  for  which 
A7  =  2  is  used.  When  A77  =  2,  that  is,  when  the  grid  is  homogeneous,  the  stability 
region  is  simply  given  as  a  linear  line  shown  as  the  blue  straight  line  in  the  figure. 
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As  the  domain  size  ratio  between  Q1  and  Qn  increases  the  stability  region  becomes 
broader.  The  green  and  red  dotted  lines  in  the  figure  represent  the  stability  region  for 
A11  =  4  and  A11  =  10,  respectively.  The  black  solid  line  represents  the  limit  of  the 
stability  region,  i.e.  for  A11  — >  oo.  The  limit  line  is  given  by  < 

0  and  is  independent  of  the  value  of  A11 . 


If  T\  —  T4  —  0,  the  penalty  interface  conditions  are  basically  the  same  as  the  upwind 
methods.  If  T\  =  r2  and  r3  =  r4  then  the  scheme  does  not  split  the  flux  into  the 
positive  and  negative  ones  but  uses  the  flux  itself  in  the  penalty  terms. 


It  is  important  to  consider  the  averaging  method  with  the  inhomogeneous  grid  system 
since  it  is  the  popular  and  simplest  interface  conditions.  For  the  averaging  method,  the 
continuity  of  q  at  the  interface  is  ensured.  The  averaging  method  is  a  special  case  of 
the  penalty  method.  The  averaging  method  has  been  effectively  used  in  the  previous 
work  [11]  for  the  homogeneous  grid  system.  Since  we  suppose  that  the  boundary 
operators  B  and  the  spectral  vanishing-viscosity  SV  ensure  the  stability  at  the  outer 
boundaries  x  =  x^,  and  Xr,  we  consider  only  the  contributions  from  the  interface  at 
x  =  0. 


Now  we  consider  the  following  penalty  scheme  for  the  averaging  method 


-g 2f  +  ~ ^  -  g(tf(o,t))]  + 

r2QN(x)[f-(qIN( 0,t))  -  /~(q£(0,t))], 

dj§  +  dV^"')  =  M[/,+«(o,  t))  -  /„+(«U0,  *))]  + 

uQM(x)[f-(q^(0,t))  -  /“(?;( 0,t))],  (13) 


where  denotes  the  derivative  with  respective  to  x.  The  averaging  method 
considered  here  is  for  the  case  that  each  subdomain  can  have  different  polynomial 
orders  such  as  N  and  M  for  Q1  and  Qn ,  respectively.  In  [11],  N  =  M  has  been  used 
for  the  averaging. 


[Averaging,  [11]]  The  scheme  (13)  is  the  averaging  method  if 


n  =  r2  =  r3  =  r4  = 


(14) 
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Proof:  If  (14)  is  satisfied  then  (13)  becomes,  at  x  =  0 


QQn  =  dq^_ 
dt  dt 


-jjj.  (|  (4/te)  +  -d 


(15) 


If  (14)  is  satisfied,  the  scheme  (13)  is  conservative  for  any  N  and  M  and  A1  and 
A11 .  Proof:  By  multiplying  (13)  by  the  weight  vectors  u>rN  and  tu1^  and  using  the 
conditions  (14)  we  have 


r°  d q:N 


l  XL 


dt 


dx  + 

Jo 


xR  Qq 


dt 


dx  = 


dft 


'XL 


dx 


dx 


fxR  df 


dx 


dx 


A1  A11 

+  |[/,te(0,t))T<  -  fx{q“( 0,f))— <] 


A11 

iifxWui  o,t))—^o 


A1 


fxWN{  o,t)v;;. 


Thus  ignoring  the  outer  boundaries,  the  RHS  of  the  above  equation  becomes 


RHS  =  -/'(<&  (0,*))  +  /"(?£(  0,t))  =  0.  (16) 

Here  we  used  the  fact  that  qIN( 0,t)  =  g"(0,t)  and  f^(qIN( 0,f))  =  /”(9m(0, t))  from 
Theorem  2.1.  □ 


The  conditions  of  the  penalty  parameters  obtained  above  are  independent  of  the 
orders  of  each  subdomain.  Moreover,  they  are  independent  of  the  domain  size  as 

well. 


With  (14)  the  scheme  (13)  is  not  necessarily  stable  in  general. 


Proof:  By  taking  the  discrete  norms,  the  energy  of  (13)  E{t)  =  ^7  q2(x,t)dx  + 
-^77  I0XR  q2(x,t)dx  becomes  without  the  boundary  terms 


1  dE(t ) 

2  dt 


+  (?w(M)  -  Qm( M))  [Aqi  -  Aq'J 
2(A7  -  A11) 


A1  A11 


(qN(0,t))  A  qN(0,t). 
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where  we  use  that  qrN(0,  t)  =  q^(Q,t)  and  r\  =  r2  =  r3  =  r4  =  7.  Thus  if  A7  7^  A77, 
the  RHS  does  not  necessarily  non-positive.  The  RHS,  however  vanishes  if  A7  =  A77 
so  that  E(t)  —  E( 0).  □ 


The  continuity  of  q  at  the  interface  is  ensured  by  using  the  averaging  method  (13). 
This,  however,  does  not  necessarily  imply  that  the  first  derivative  is  also  continuous 
at  the  interface.  In  general  the  first  derivative  is  discontinuous. 


2.2  Weighted  spectral  penalty  method 


In  the  previous  section,  it  has  been  shown  that  the  conservative  and  stable  penalty 
method  can  be  constructed  for  fully  inhomogeneous  grids,  i.e.  N  7^  M  and  A7  7^  A77 
except  the  averaging  method.  The  conditions  obtained  in  the  previous  section  are 
only  necessary  conditions.  For  example,  the  stability  conditions,  (7)  and  (8)  say  that 
the  scheme  (3)  is  stable  if  T\  =  r2  =  and  73  =  74  =  —  yyr-  These  conditions,  albeit 
stable,  can  yield  the  nonphysical  reflecting  solutions  at  the  inhomogeneous  domain 
interfaces  because  the  positive  and  negative  fluxes  are  equally  penalized  as  we  will 
show  in  this  section.  In  this  section,  the  weighted  spectral  penalty  method  for  the 
inhomogeneous  grid  is  introduced  to  reduce  the  nonphysical  modes  at  the  domain 
interfaces. 


With  the  weighted  spectral  penalty  method,  the  incoming  or  outgoing  characteristics 
are  penalized  with  different  weights  if  the  inhomogeneous  domain  system  is  considered 
such  that  the  incoming  fluxes  are  penalized  with  the  larger  values  of  the  penalty 
parameters  than  the  outgoing  fluxes.  In  the  Legendre  spectral  penalty  equation  (3), 
the  weighted  spectral  penalty  method  for  f V  exploits 

N  >  M, 


and  for  Q11 


r3\  >  t4 


The  numerical  simulation  results  of  supersonic  reactive  cavity  flow  presented  in  this 
paper  show  that  the  upwind  characteristic  interface  conditions  are  not  enough  to  en¬ 
sure  the  smooth  solutions  across  the  interfaces.  We  will  show  in  the  following  sections 
that  by  weighting  the  incoming  fluxes  against  the  outgoing  fluxes,  the  nonphysical 
modes  at  the  domain  interfaces  can  be  reduced.  The  weight,  however,  can  not  be 
arbitrarily  large  due  to  the  CFL  restriction.  In  practice,  we  use  the  fixed  weight  for 
each  penalty  parameter.  Since  the  problem  considered  in  this  paper  is  highly  non- 
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linear,  the  fixed  weight  for  Vt  >  0  may  not  be  enough  to  prevent  the  growth  at  the 
interfaces.  We  use  the  local  spectral  vanishing  viscosity  method  with  the  weighted 
penalty  method  to  prevent  any  growth  at  the  interfaces. 


2.2.1  Reflection  coefficients  of  the  weighted  penalty  interface  conditions 


In  order  to  explain  how  the  weighted  spectral  penalty  method  can  reduce  the  non¬ 
physical  reflection  modes,  the  reflection  coefficients  analysis  is  used.  Here  the  reflection 
modes  are  obtained  due  to  the  different  resolution  of  each  subdomain  across  the 
domain  interfaces  although  the  solution  can  be  continuous.  Consider  the  following 
simple  linear  hyperbolic  equation 


qt  +  (. Fq)x  =  0,  q  :  M  x 


x  G 


t  >  0, 


(17) 


where 


Q  = 


H 

,  F  = 

/ 0  l\ 

w 

yioj 

(18) 


The  same  equation  has  been  considered  to  show  the  reflecting  modes  at  the  domain 
interfaces  with  the  spectral  Galerkin  method  in  [15].  We  seek  a  wave  solution  such 
that 

q(x,  t)  —  exp(iut)q(x),  x  &  [—2,2],  t>  0.  (19) 

Plugging  the  wave  solution  into  (18)  yields 


q(x)  =  Aqi  exp(— iojx)  +  Bq2exp(ioux), 
where  q\  =  (1, 1)T  and  q2  =  (1,  —  1)T. 


(20) 


Suppose  that  we  have  two  subdomains  fP  =  [—2,0],  and  ft11  =  [0,2],  Let  B±  be  the 
boundary  operators  at  the  end  points,  i.e.  x  =  —  2  and  x  =  2.  And  f+  and  f~  are 
f±  =  Ffq  =  SA±S~1q  ; 


f+  = 


u  +  v 

U  +  V 


r  = 


(21) 


Here  we  assume  that  the  boundary  operator  B  is  taken  properly  such  that  this 
treatment  does  not  destroy  the  global  stability  and  there  is  no  reflection  from  the 
boundaries.  In  other  words,  we  assume  that  we  have  the  perfect  and  stable  absorbing 
boundary  operator  at  x  —  ±2.  Plugging  the  wave  solutions  into  the  Legendre  spectral 
method  (3),  we  have  the  following  linear  system  at  the  interface,  i.e.  at  x  —  0, 
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niA1  -  A11)  -  t2[B!  -  B11)  =  0,  r3(An  -  A1)  -  T4(Bn  -  B1)  =  0, 
niA1  -  A11)  +  t2(B!  -  B11)  =  0,  t3{Au  -  A1)  +  u(Bn  -  B1)  =  0. 

The  above  linear  systems  can  be  rewritten  in  the  matrix  formWX  =  ZX  with 


(  Ai  > 

B1 

A11 


\BU  ) 


W 


t  —Ti  Ti  0  0  ^ 

— r2  r2  0  0 

0  0  —t3  t3 

\  0  0  — r4  r4  j 


^00  —  Ti  Tx  ^ 

0  0  -r2  r2 

-r3  r3  0  0 

^  -r4  r4  0  0  y 


The  system  is  not  well-posed  as  det(W  —  Z)  =  0.  In  fact  this  linear  system  can  be 
solved  by  taking  into  account  that  A11  and  B11  are  considered  as  the  given  boundary 
values  for  the  solution  of  VL1  and  A1  and  B1  of  Q,11  in  the  real  computation.  In  order  to 
look  at  the  possible  reflecting  modes  at  the  interface,  assume  that  the  non-zero  wave 
solutions  are  locally  defined  only  in  IT,  that  is,  the  solution  is  compactly  supported 
in  Q1  such  that  it  vanishes  in  Thus  A11  =  B11  =  0.  We  take  the  ratio  of  the 
coefficients  at  the  interface  using 


TiA1  —  t2B!  =  0,  TiA1  +  t2B!  =  0. 


Define  the  reflection  coefficients  Rq  at  x  —  0  for  Since  A  is  corresponding  to 
the  outgoing  flux  and  B  to  the  incoming  flux  at  x  =  0,  respectively,  the  reflection 
coefficients  R0  is  given  by 


Rq  — 


B1 

Tl 

A1 

T2 

(22) 


Table  I 

Interface  conditions  and  reflection  coefficients. 


Case 

Interface  conditions 

Reflection  coefficients 

T-l 

Tl  =  0 

o 

II 

o 

T-2 

T\  =  T2  =  T 

Rq  =  1 

T-3 

T\  <  T2 

Rq  1 

T-4 

Tl  »  T2 

Rq  1 

We  shall  consider  four  different  cases  as  given  in  table  I.  By  definition,  there  is  no 
reflection  at  the  interface  for  the  case  T-l,  the  upwind  method.  For  the  case  T-2,  we 
do  not  split  the  flux 

Bf+  +  r2/“  =  r(/+  +  /")  =  rf.  (23) 
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Consequently,  the  reflection  is  obvious  although  the  method  is  simple.  The  case  T-3 
is  the  weighted  penalty  method.  The  case  T-4  weights  the  outgoing  flux  against  the 
incoming  flux. 


For  illustration,  consider  a  step  function,  i.e. 


liefl' 

U  =  V  =  < 

o  x  g  n11 


Figure  3  shows  the  solutions  after  one  time  integration  at  t  =  At  =  0.0001  for  each 
case  in  Table  1.  Ti  =  r2  =  are  used  for  the  case  T-2,  r2  ~  0(N3)  and  T\  ~  0(N2) 
are  used  for  the  case  T-3  and  T\  ~  0(N3)  and  r2  ~  0(N 2)  are  used  for  the  case  T-4. 
Note  the  different  behaviors  of  the  case  T-3  and  the  case  T-4.  For  the  case  T-3  the 
penalty  parameters  associated  with  the  incoming  flux  are  weighted  while  the  outgoing 
flux  are  weighted  for  the  case  T-4.  The  solution  for  the  case  T-3  does  not  show  the 
overshoot  at  the  interface  for  Q1 .  For  the  case  T-4,  the  solution  of  fF  at  interface 
shows  the  overshoot.  The  behaviors  of  the  interface  solution  of  Q11  for  each  case  can 
be  explained  by  taking  into  account  that  the  scheme  is  in  fact  conservative. 


(a)  T-l:  T\  =  0 


(b)  T-2:  n  =  r2  =  zj- 


(c)  T-3:  t2  ~  0(N3)  (d)  T-4:  n  ~  0(N3) 


Fig.  3.  Shock  calculation  at  the  first  time  step 
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2.2.2  Reflection  and  instability 


We  consider  more  numerical  examples  to  confirm  the  performance  of  the  weighted 
spectral  penalty  method. 


Homogeneous  grid:  First  consider  the  same  wave  equation  with  the  boundary  condi¬ 
tions 


q  =  (u,v)T,  f  =  Fq  =  (v,  u)T, 
u(x,  t )  —  v{x ,  t)  =  0,  x  =  4 

Bq\=l 

I  u(x,  t )  +  v(x,  t)  =  0,  x  =  0 

where  B  is  the  boundary  operator.  Here  we  again  consider  two  subdomains  VI1  =  [0,  2] 
and  fW  =  [2,4]  with  the  same  polynomial  order  N.  For  the  boundary  conditions  at 
x  =  0  and  x  =  4,  we  use  the  non- reflecting  boundary  conditions.  The  spectral  penalty 
method  (3)  described  in  the  above  sections  is  used  with  the  Legendre  polynomials. 


(24) 

(25) 


We  denote  M-A,  M-UW,  M-NFS  and  M-WP  by  the  averaging,  upwind,  no-flux  split¬ 
ting  and  the  weighted  penalty  methods,  respectively.  They  are  listed  in  the  following 
table  below. 


Method 

Interface  Conditions 

Remark 

M-A 

T\  =  r2  =  r3  =  r4  =  \ 

Averaging  Method 

M-UW 

n  =  74  =  0,72  =  -73  =  ^ 

Upwind  Method 

M-NFS 

Tl  =  T2  =  -r3  =  -74  =  2^ 

No-flux-Splitting  Method 

M-WP 

0(7-2)  =  0(73)  ~  N3,  ri/0/r4 

Weighted  Penalty  Method 

Note  that  all  four  methods  satisfy  the  stability  condition  derived  in  (8)  as  each  sub- 
domain  has  the  same  polynomial  order  N  and  domain  length  A. 


The  CFL  condition  is  given  by 


mm 


At 

Axi 


<  C 


(26) 
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where  Cfl  is  a  positive  constant  and  taken  to  make  At  is  small  enough  for  every 
case. 

Table  II 

The  maximum  error  (L0 c)  for  M-A,  M-UW,  M-NFS  and  M-WP  at  time  t  =  1.5. 


N\Method 

M-A 

M-UW 

M-NFS 

M-WP 

4 

0.67 

0.82 

0.13E+01 

0.74 

8 

0.65E-02 

0.28E-02 

0.82E-02 

0.57E-02 

16 

0.34E-09 

0.38E-09 

0.65E-09 

0.29E-09 

32 

0.52E-11 

0.52E-11 

0.52E-11 

0.52E-11 

64 

0.54E-11 

0.51E-11 

0.54E-11 

0.50E-11 

Table  II  shows  the  L ^  error  for  each  method.  The  overall  performance  is  almost 
the  same  for  each  method  while  the  M-WP  performs  slightly  better  than  the  other 
methods. 


Inhomogeneous  grid 


Now  we  consider  the  same  problem  with  3  subdomains,  each  of  them  having  the 
same  domain  intervals,  i.e.  fl1  =  [0,  2] ,  O2  =  [2,4]  and  fl3  =  [4,6].  For  these  three 
subdomains,  consider  the  following  two  different  cases: 


Case 

Grid  resolution 

C-l  (homogeneous) 

C-2  (inhomogeneous) 

N\  =  IV2  =  A3  =  8 

Ni  =  8,  A2  =  32,  A3  =  8 

Note  that  every  method  satisfies  the  stability  condition  for  the  case  C-l.  For  the 
case  C-2,  only  the  M-A,  which  is  the  averaging  method,  does  not  satisfy  the  stabil¬ 
ity  condition  as  explained  in  Section  2.1,  but  the  M-NFS  still  satisfies  the  stability 
condition. 


Table  III  shows  the  L ^  error  for  each  case  for  the  whole  domain.  As  shown  in  the 
table,  the  M-A  and  the  M-NFS  show  the  instability  at  the  interfaces  for  the  case  C-2. 
Table  IV  shows  the  L <*,  error  for  hi2  for  each  case.  The  L <*,  error  for  fl2  is  less  than 
that  of  global  error  for  the  M-UW  and  the  M-WP  of  the  case  C-2  because  the 
higher  polynomial  order  of  N  is  used. 
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Table  III 

The  maximum  error  (L0 c)  for  the  case  C-l  and  the  case  C-2  at  time  t  =  1.5. 


Case\Method 

M-A 

M-UW 

M-NFS 

M-WP 

C-l 

0.83E-02 

0.29E-02 

0.11E-01 

0.52E-02 

C-2 

unstable 

0.25E-02 

unstable 

0.65E-02 

Table  IV 

The  maximum  error  (Loo)  °f  fl2  for  C-l  and  C-2  at  time  t  =  1.5. 


Case\Method 

M-A 

M-UW 

M-NFS 

M-WP 

C-l 

0.83E-02 

0.29E-02 

0.10E-01 

0.51E-02 

C-2 

unstable 

0.13E-02 

unstable 

0.51E-02 

The  instabilities  of  the  M-A  and  the  M-NFS  are  illustrated  in  Fig.  4.  In  the  figure,  the 
top  figures  represent  the  M-A  at  t  =  0.4  and  the  bottom  figures  the  M-NFS  at  t  —  1.5. 
The  left  figures  show  u+v  and  the  right  u  —  v.  As  shown  in  the  figures,  the  locations  of 
the  instability  are  different  for  the  M-A  and  the  M-NFS.  Since  N2  >  Ni  =  N3,  there 
exist  higher  modes  in  f I2  which  do  not  appear  in  the  approximations  for  both  fl1  and 
hi3.  For  the  M-A,  i.e.  the  averaging  method,  the  figure  indicates  that  the  instability 
occurs  at  the  interface  of  fl1  and  hi2  for  u  +  v,  and  the  instability  at  the  interface  of 
hi2  and  hi3  for  u  —  v.  This  implies  that  the  sudden  growth  at  the  interface  occurs  when 
the  characteristic  of  the  lower  modes  enters  the  subdomain  where  the  higher  modes 
appear  in  the  approximation.  For  the  M-NFS,  the  growth  occurs  at  different  locations. 
If  the  outgoing  characteristic  is  approximated  with  the  higher  modes,  such  modes  in 
the  approximation  are  reflected  as  if  the  adjacent  subdomain  plays  a  role  as  a  wall 
boundary.  The  subdomain  fl2  yields  a  free  boundary  condition  for  the  lower  mode 
wave  solutions  entering  hi2 .  No  significant  growth  at  the  interface  is  not  observed  for 
the  weighted  penalty  method. 


2.2.3  Adaptive  super-viscosity  method 


As  mentioned  above,  we  use  the  fixed  weight  for  each  penalty  parameter  for  the  real 
computation  and  this  is  not  enough  for  preventing  the  growth  at  the  interface  if 
the  grid  system  is  highly  inhomogeneous.  The  super- viscosity  method  is  used  at  the 
interface  when  the  weighted  penalty  method  fails  to  reduce  the  nonphysical  growing 
modes.  The  super-viscosity  method  is  only  applied  when  and  where  it  is  necessary. 
Since  the  nonphysical  reflecting  or  growing  modes  are  mainly  contained  in  the  higher 
modes,  the  super-viscosity  can  be  used  to  remove  such  higher  modes  locally  at  the 
interface.  This  technique  is  critical  for  the  stabilization  of  the  multi-domain  spectral 
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(a)  M-A,  (u+v) 


(b)  M-A,  (u-v) 


(c)  M-NFS,  (u+v)  (d)  M-NFS,  (u-v) 


: 

6 

- 

2 

1  0 

-2 

: 

-6 

- 

Fig.  4.  Solutions  for  the  M-A  at  time  t  =  0.4  and  the  M-NFS  at  time  t  =  1.5.  The  u  +  v 
and  u  —  v  characteristic  waves  are  shown  in  the  left  and  right  columns,  respectively. 

penalty  method  for  the  nonlinear  problems. 

For  the  numerical  experiment,  we  revisit  the  simple  wave  problem  used  in  [5], 


du  |  du  n 

dt  ‘r  dx  ~  U 

x  E  Cl  =  (a,b),  t  >  0 

u(a,t)  =  uL(t ) 

t  >  0 

(27) 

u(x,  0)  =  COs(7Tx) 

x  E 

We  seek  an  approximation  of  u  with  two  subdomains,  fl1  =  [0,  2]  and  fl2  =  [2, 4]  with 

■V,  /  ,\2. 

Two  cases  of  grid  resolution  N  in  each  subdomain  are  examined  for  the  M-A  and  the 
M-NFS  : 


Case 

Grid  resolution 

C-3 

C-4 

iVi  =  32,  N2  =  8 

Ni  =  8  ,N2  =  32 
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Table  V  shows  the  errors  of  the  case  C-3  and  the  C-4.  The  table  shows  that  the 
methods  are  stable  in  either  case  except  the  M-A  and  the  M-NFS.  Both  of  them  are 
unstable  if  N\  <  TV2  and  TV2  <  TVi,  respectively. 

Table  V 


maximum  error 

[Loo)  for  C-3  and  C-4  at  time  t  = 

1.5. 

Case/Method 

M-A 

M-UW 

M-NFS 

M-WP 

C-3 

0.13E-01 

0.46E-02 

4.89  (Unstable) 

0.49E-02 

C-4 

0.54E+17  (Unstable) 

0.46E-02 

0.461E-02 

0.46E-02 

To  avoid  the  growth  in  time  the  artificial  super  vanishing  viscosity(SV)  term  is  added 
at  the  interface: 


duN 

~dt 


+ 


duN 

dx 


1 

N2s~i 


d_ 

dx 


(1-x2) 


un  +  PT 


x  e 


(28) 


where  s  is  a  positive  integer  growing  with  TV  [16,20]  and  PT  denotes  the  penalty 
term.  This  SV  method  is  equivalent  to  the  filtering  method  [16,20].  The  exponen¬ 
tial  filter  method  with  the  filtering  order  7  is  used  for  the  numerical  experiment. 
The  exponential  filter  function  <r(k)  and  the  filtering  order  7  are  defined  as  cr(k)  = 
exp (—€M(k/N)'y)  where  k  is  the  mode  number  k  —  0,  •  •  •  ,  TV  and  TV  is  the  polynomial 
order  used.  The  positive  constant  6m  is  chosen  such  that  <r(TV)  becomes  machine  zero. 
Typically  ~  32. 


Tabic  VI  shows  the  results  for  the  M-A  and  the  M-NFS  with  the  different  orders 
of  filtering  7.  As  the  table  indicates  there  no  significant  growth  has  been  observed. 
Comparing  with  the  results  in  Table  II,  however,  we  notice  that  the  filtering  method 
permits  a  lose  of  accuracy.  For  example,  the  error  of  M-NFS/C-3(TV  =  (32,8)) 
is  0.11  x  10-1  and  M-NFS/C-4  (TV  =  (8,32))  is  0.46  x  10~2.  Table  VI  also  indicates 
that  the  methods  are  stable  even  though  7  — »  00,  which  implies  that  stability  can  be 
achieved  even  with  the  small  viscosity  added. 


Fixing  the  filtering  order  7  =  16,  Table  VII  shows  the  full  recovery  of  the  accuracy 
as  TV  increases. 


Here  we  note  that  the  weighted  penalty  method  is  proposed  for  the  ID.  For  the  2D 
problem,  one  needs  to  find  the  proper  conditions  for  the  corner  of  each  subdomain. 
Such  conditions  will  be  investigated  in  our  future  work.  The  2D  numerical  experiments 
indicate,  however,  that  the  proposed  method  with  the  SV  method  applied  at  the 
interface  and  corner  yields  a  stable  and  accurate  result  as  shown  in  the  next  sections. 
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Table  VI 

The  maximum  error  (L0 c)  for  the  reconstruction  method. 


Table  VII 

The  maximum  error  (L^)  for  the  reconstruction  method  ;  7  =  16. 


(M-A/C-4)  (U1) 

(M-NFS/C-3)  (U1) 

N 

=  4 

1.24 

1.00 

N 

=  8 

0.76E-01 

0.11E-01 

N 

=  16 

0.47E-06 

0.12E-08 

N 

=  32 

0.74E-11 

0.52E-11 

3  Injector-cavity  scramjet  system 


In  this  section,  the  proposed  weighted  spectral  penalty  method  is  applied  for  the 
approximation  of  the  supersonic  flow  interactions  for  the  cavity-injector  system  with 
the  inhomogeneous  grid  system.  To  refine  the  localized  injector  field,  we  use  the 
narrow  subdomain  for  the  injector.  The  polynomial  order  of  each  domain  is  the  same 
for  both  x—  and  y—  directions.  Thus  the  grid  inhomogeneity  in  this  case  comes  from 
the  different  domain  length. 


3.1  The  cavity  system  and  the  governing  equations 


Cavity  has  been  actively  used  as  a  flame-holder  in  scramjet  engine  (see  the  review  by 
Ben-Yakar  and  Hanson  [4]).  The  injector-cavity  system  is  illustrated  in  Figure  1.  The 
cavity  system  is  categorized  into  4  different  types  such  as  open,  closed,  transitional- 
closed  and  transitional-open  depending  on  the  length  scale  of  cavity  [4],  The  cavity 
system  with  the  length-to-depth  ratio  L/D  <  7  ~  10  is  called  an  open  cavity  as  the 
upper  shear  layer  reattaches  itself  at  the  back  face.  Under  the  shear  layer  formed 
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over  cavity,  the  flows  with  the  hydrogen  fuel  are  possibly  captured  inside  cavity  and 
generate  the  recirculation  zone.  The  generated  recirculation  interacts  with  the  shear 
layer  and  the  acoustic  waves  inside  cavity.  The  radicals  from  the  chemical  reaction 
between  the  hydrogen  and  oxygen  gases  reside  inside  cavity  and  trigger  the  auto¬ 
ignition  of  the  supersonic  engine.  In  principle,  the  more  stable  and  longer  recirculation 
is  maintained,  the  more  efficient  fuel  performance  can  be  achieved. 


The  major  question  of  the  cavity  flame-holder  system  that  needs  to  be  investigated 
is:  How  does  the  fuel  injection  interact  with  cavity  flows?  There  have  been  many 
numerical  studies  on  the  recirculation  and  stabilizations  of  the  flow  inside  cavity 
but  rarely  on  how  the  continuous  supply  of  the  fuel  can  affect  the  flow  dynamics 
inside  cavity  [2,3,7,21,22,25,26].  Since  the  injection  of  the  fuel  in  the  combustor  is 
necessary,  however,  the  injection  emerges  as  another  important  key  parameter  for  the 
optimal  configuration  of  the  cavity  flame-holders.  Both  comprehensive  laboratory  and 
numerical  experiments  have  to  be  carried  out  to  answer  the  question.  In  this  work,  we 
use  the  length-to-depth  ratio  L/D  =  4cm/lcm  =  4,  that  is,  we  use  the  open  cavity. 


The  governing  equations  are  the  compressible  2D  reactive  Navier-Stokes  equations 
with  the  chemical  source  terms  given  by 


dq  dF_  dG__dF\^  dGu 
dt  dx  dy  dx  dy 


(29) 


where  q  =  (p,  pu,  pv,  E,  pf)T  is  the  state  vector,  F  =  ( pu ,  pu2  +  P ,  puv,  ( E  +  P)u ,  pfu)T 
and  G  =  (pv,  puv ,  pv2  +  P,(E  +  P)v,  piv)T  the  inviscid  fluxes,  Fu  and  Gu  the  viscous 
fluxes  and  C  the  chemical  source  term,  respectively.  Here  p,u,v,  E,  P,  and  f  denote 
the  density,  the  velocity  in  ^-direction  and  the  velocity  in  ^/-direction,  the  total  energy, 
the  pressure  and  the  mass  fraction  vector,  respectively.  The  chemical  model  uses  four 
chemical  species,  H2,  02,^0  and  N2  with  the  reversible  chemical  reaction  between 
hydrogen  and  oxygen  gases  given  by 


2H2  +  02  ^  2H20.  (30) 

A  modified  Arrhenius  law  gives  the  equilibrium  reaction  rate  ke,  the  forward  reaction 
rate  kf  and  the  backward  reaction  rate  kb  as 


ke  =  AeT exp(4.60517(i7e/T  -  2.915)), 
kf  =  Af  exp (-Ef/(RT))} 
k\)  kf  j  ke  7 

where  Ee  =  12925,  and  Ef  =  7200  are  the  activation  energy  and  Ae  =  83.006156,  and 
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Af  =  5.541  x  1014  are  the  frequency  factors.  R  is  the  universal  gas  constant.  Each 
chemical  species  has  different  dynamical  viscosity  /q  based  on  the  Sutherland’s  law. 
The  mixture  viscosity  fi  is  obtained  according  to  the  Wilke’s  law  [24],  The  Prandtl 
number  Pr  and  the  Schmidt  number  Sc  are  taken  to  be  0.72  and  0.22  respectively 
for  the  normal  air.  The  equation  of  state  is  given  by  the  assumption  of  the  perfect 
gas  law.  Detailed  formulation  of  the  equations  can  be  found  in  [9,11] . 


With  L  fixed  we  consider  four  different  angles  of  the  aft  wall,  i.e.  90°  and  30°.  For  the 
fluid  conditions,  the  free  stream  Mach  number  M  =  1.91,  total  pressure  P  =  2.82atm, 
total  temperature  T  =  830. 6K  and  normalized  Reynolds  number  Re  =  3.9  x  107m  1 . 
Note  that  the  Reynolds  number  is  normalized  and  has  the  unit  of  1/ [length],  and 
that  the  Reynolds  number  based  on  the  cavity  dimensions  is  about  O(105).  The 
boundary  layer  thickness  scale  is  S  =  5  x  10_4m,  and  the  wall  temperature  is  Tw  = 
460.7835K.  For  more  detailed  physical  configuration  and  its  explanation,  we  refer 
[11],  The  hydrogen  fuel  is  injected  1.5cm  ahead  of  the  cavity  with  the  injection  Mach 
number  M  =  1  (see  Figure  1).  The  numerical  experiments  are  conducted  with  two 
different  sizes  of  the  injector  diameter,  d  =  2mm  and  d  =  2cm  to  investigate  the 
effect  of  the  injector- channel  flow  interactions  on  the  development  of  the  shear  layer 
over  cavity.  The  fuel  is  injected  into  the  channel  flow  with  the  direction  normal  to 
the  base  wall.  The  total  pressure  and  the  total  temperature  of  the  hydrogen  jet  are 
2.828522atm  and  830. 6K  respectively. 


3.2  Grid  inhomogeneity  and  the  weighted  penalty  method 


To  deal  with  the  grid  inhomogeneity  we  use  the  weighted  penalty  method  described 
in  Section  3. 


The  weighted  penalty  method  is  based  on  the  characteristic  decomposition  and  it 
does  not  modify  the  stability  conditions  associated  with  Au  ■  q  and  Au  ■  dq  for  the 
Navier-Stokes  equations  in  [11], 


Figure  5  shows  the  effect  of  the  grid  inhomogeneity  on  the  solution.  The  subdomain 
containing  the  injector  has  a  smaller  domain  length  than  the  channel  subdomains. 
The  left  figure  shows  the  solution  based  on  the  averaging  method  (M-A)  and  the  right 
shows  the  solution  based  on  the  weighted  penalty  method  (M-WP).  The  figures  clearly 
show  that  the  averaging  interface  condition  (M-A)  yields  a  nonphysical  concentration 
near  the  domain  interface  while  the  weighted  penalty  interface  condition  (M-WP) 
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yields  smooth  solutions  across  the  interfaces.  Figure  6  shows  the  results  near  the 
injector  subdomain  with  the  weighted  penalty  method.  The  figure  shows  the  flow 
streamline  near  the  injector.  The  figures  show  that  there  is  no  significant  reflection 
at  the  injector  subdomain  interfaces.  It  is  also  shown  that  the  flow  fields  are  well 
resolved  with  the  weighted  penalty  method.  The  small  recirculation  formed  in  front 
of  the  injector  is  clearly  seen.  Such  recirculation  is  physically  formed  due  to  the 
interaction  between  the  incoming  channel  flow  and  the  hydrogen  jet  with  the  no-slip 
boundary  condition  at  the  wall  [3,4], 


Fig.  5.  Density  contour  of  injector-cavity-channel  flow  by  reactive  Navier-Stokes  equations 
for  the  normal  cavity  flame-holder.  Each  subdomain  has  the  same  polynomial  order  for  the 
approximation,  that  is,  each  subdomain  has  N  grid  points  both  in  x—  and  y— directions 
but  the  different  subdomain  length.  The  left  figure  shows  the  solution  using  the  averaging 
interface  conditions  (M-A)  and  the  right  figure  shows  the  solution  using  the  weighted  penalty 
interface  conditions  (M-WP).  The  weighted  penalty  interface  condition  method  considerably 
reduces  the  nonphysical  density  concentration  near  the  interfaces  seen  in  the  left  figure. 


Fig.  6.  The  recirculation  zones  formed  in  front  of  the  hydrogen  jet:  the  flow  streamlines  are 
given  with  the  hydrogen  jet  contour  for  the  narrow  injector-cavity  system  at  t  =  0.225ms. 
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3.3  Shear  layer  interactions 


One  of  the  major  effects  on  the  stability  of  the  recirculation  zone  is  by  the  shear 
layer  over  cavity.  In  [4]  (also  see  references  therein),  several  effects  on  the  shear  layer 
formation  and  its  interaction  with  the  cavity  have  been  discussed  including  the  lo¬ 
cation,  size  and  the  total  number  of  injectors.  For  the  numerical  experiments,  we 
consider  two  different  injectors,  the  narrow  and  broad  injectors.  Figure  7  shows  the 
water  contours  for  both  cases.  By  placing  the  injector  ahead  of  the  cavity  front  wall, 
the  pressure  fluctuations  are  reduced  and  the  sharp  gradients  found  near  the  corner 
of  the  aft  wall  are  also  weakened  as  the  shear  layer  is  being  developed.  The  figures 
show  that  the  broader  injector  has  more  enhanced  shear  layer  growth  over  the  cavity 
than  the  narrow  jet.  However,  the  pressure  profiles  in  Figures  9  and  ??  indicate  that 
the  pressure  oscillations  can  be  more  attenuated  with  the  narrow  injector  than  the 
broad  injector.  This  implies  that  there  exists  an  optimal  size  of  the  injector  with  the 
fixed  location  of  the  injector  from  the  front  cavity  wall  that  minimizes  the  pressure 
fluctuations  and  maximizes  the  stability  of  the  recirculation  zones  inside  cavity.  Both 
the  broad  and  narrow  injector  systems  also  show  that  they  have  weaker  flow  gradi¬ 
ents  near  the  aft  wall  than  the  flow  gradients  obtained  in  our  previous  work  without 
the  injector.  The  injection  angle  is  normal  to  the  wall  but  different  injection  angles 
can  be  used.  In  [3,4],  it  has  been  discussed  that  the  angled  injector  such  as  30°  or 
60°  can  further  weaken  the  possible  bow  shock  found  at  the  aft  wall.  Figure  8  shows 
some  detailed  differences  of  the  water  and  hydrogen  profiles  between  the  normal  and 
recessed  cavities  at  t  =  3.48ms  for  the  narrow  injector  system. 


3-4  Pressure  fluctuations 


In  [11],  we  considered  the  cold  and  reactive  flows  without  the  hydrogen  injector  and 
showed  that  the  pressure  fluctuations  inside  cavity  can  be  considerably  reduced  if  the 
aft  wall  is  slanted.  Consequently  this  helps  more  stable  recirculation  inside  cavity  to 
be  developed.  The  generated  acoustic  waves  disturbing  the  recirculation  are  reflected 
back  to  the  shear  layer  due  to  the  slantness  of  the  rear  wall.  Similar  results  are 
found  in  the  cavity  system  with  the  injection  fields.  Figure  9  shows  the  pressure 
fluctuation  history  profiles  for  the  normal  (90°,  top  figures)  and  slanted  wall  (30°, 
bottom  figures)  cases  for  the  broad  (left)  and  narrow  (right)  injectors.  In  the  figures, 
the  pressure  fluctuations  are  measured  to  tU0/D  ~  150  but  plotted  in  the  same 
scale  used  in  Figure  6  of  [11]  for  the  comparison.  The  pressures  are  measured  at 
the  center  of  cavity.  For  the  broad  injector-cavity  system,  it  is  clearly  shown  that 
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Injector 


Injector 


Fig.  7.  Water  production  and  transport  for  90°  and  30°  cavity  walls:  the  upper  figures  shows 
the  broad  injector-cavity  system  and  the  bottom  figures  the  water  contours  of  the  narrow 
injector-cavity  system  at  t  =  0.225ms,  respectively.  There  are  50  contour  levels  ranged  from 
-0.001  to  0.23. 


Fig.  8.  Water  and  hydrogen  density  profiles  at  t  =  3. 48ms  for  the  narrow  injector-cavity 
system  with  30°  and  90°  aft  walls.  The  top  figures  show  the  water  density  contours  and  the 
bottom  figures  the  hydrogen  density  contours. 

the  pressure  fluctuations  are  much  attenuated  for  the  lower  aft  wall  case  and  these 
features  are  similar  to  those  for  the  non-reactive  cold  flow  cases.  For  the  narrow 
injector-cavity  system,  the  pressure  fluctuations  for  both  90°  and  30°  wall  cavities 
are  highly  attenuated  compared  to  those  for  the  broad  injector-cavity  system.  The 
differences  of  the  pressure  fluctuations  between  30°  and  90°  are  not  significant,  but 
one  can  observe  that  the  lower  angled  wall  cavity  has  less  pressure  fluctuations  than 
the  normal  wall  cavity.  These  results  are  similar  to  those  for  the  reactive  flow  cases 
without  the  injection  fields.  Note  that  the  pressure  fluctuations  of  the  normal  wall 
cavity  system  are  also  much  attenuated  compared  to  the  pressure  fluctuation  of  the 
normal  wall  cavity  system  without  the  injection  fields.  The  injector  held  in  front  of 
the  cavity  increases  the  stability  of  the  recirculation  inside  cavity. 
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Fig.  9.  Pressure  history  for  broad  (left)  and  narrow  (right)  injector-cavity  flows  at  the  center 
of  cavity.  For  the  broad  jet,  the  diameter  is  d  =  2cm  and  d  =  2mm  for  the  narrow  jet.  Each 
panel  shows  the  case  of  90°  and  30°  cavity  walls  from  top  to  bottom. 

4  Summary 


In  this  research,  the  direct  numerical  simulation  of  the  supersonic  injector-cavity 
scramjet  system  has  been  carried  out  with  the  multi-domain  spectral  penalty  method 
with  the  inhomogeneous  grid  system.  In  order  to  minimize  the  development  of  the 
nonphysical  modes  generated  at  the  domain  interfaces  due  to  the  grid  inhomogeneity, 
we  first  described  the  stable  and  conservative  interface  conditions  of  the  multi-domain 
spectral  penalty  method  for  the  inhomogeneous  grid  system.  For  general  inhomoge¬ 
neous  grid  system,  it  is  shown  that  it  is  possible  to  construct  a  stable  and  conservative 
spectral  penalty  method.  For  the  averaging  method,  it  is  also  shown  that  the  conser- 
vativity  can  be  preserved  but  the  stability  is  not  maintained  in  general.  The  weighted 
penalty  interface  conditions  is  then  proposed  to  minimize  the  non-physical  effect  at 
the  inhomogeneous  grid  interfaces.  The  weighted  penalty  method  gives  more  weight 
to  the  incoming  fluxes  than  the  outgoing  fluxes.  For  the  numerical  experiments,  we 
use  the  fixed  weight  for  all  time.  The  weight,  however,  can  be  adaptively  determined 
depending  on  the  flow  conditions.  Such  adaptivity  will  be  investigated  in  our  future 
work.  The  weight  penalty  method  reduces  the  nonphysical  growth  at  the  domain  in¬ 
terface  considerably  but  not  completely.  The  adaptive  Altering  method  is  used  with 
the  weighted  penalty  method  to  stabilize  any  growth  at  the  interface.  The  adaptive 
Altering  is  only  applied  at  a  small  number  of  points  at  the  interface.  The  direct  nu¬ 
merical  simulation  shows  that  the  proposed  method  successfully  yields  the  stable  and 
accurate  approximation  of  the  injector-cavity  Aows  with  the  inhomogeneous  grids.  It 
is  qualitatively  shown  that  the  recessed  cavity  yields  a  better  performance  of  the  pres¬ 
sure  Auctuation  reduction  and  enhances  the  stability  of  the  recirculation  zones  inside 
cavity.  The  injector  located  in  front  of  the  cavity  also  reduces  the  pressure  Auctuations 
inside  cavity.  More  detailed  geometric  conAgurations  maximizing  the  attenuation  of 
the  pressure  Auctuations  and  the  stability  of  the  recirculation  inside  cavity  will  be 
investigated  in  our  future  work.  The  future  research  work  will  also  center  around  the 
development  of  the  3D  spectral  penalty  method  with  the  weighted  penalty  conditions. 
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5  Introduction 


In  this  article  we  extend  the  one- dimensional  Multi-Domain  Hybrid  Spectral- WENO 
Method  (Hybrid)  for  Hyperbolic  Conservation  Laws  [7]  to  two- dimensions  in  space 
and  apply  it  to  the  classical  Shock- Vortex  interaction  and  Richtmyer-Meshkov  insta¬ 
bilities  problem.  The  general  idea  of  the  Hybrid  method  is  to  use  a  multi-domain 
framework  in  order  to  apply  convenient  spatial  discretizations  to  the  smooth  and 
rough  parts  of  the  numerical  solution.  Shocks  and  high  gradients  are  kept  at  WENO 
subdomains,  while  complex,  but  still  smooth,  details  of  the  solution  are  treated  within 
spectral  subdomains.  Numerical  efficiency  is  increased  with  respect  to  the  classical 
spectral  and  WENO  methods:  Postprocessing  techniques  of  the  spectral  method  ap¬ 
proach  of  shocks  [15,20]  are  avoided,  since  no  Gibbs  phenomenon  will  occur,  and  the 
expensive  characteristic  decompositions  and  projections  of  the  WENO  method  are 
skipped  at  the  smooth  parts  of  the  solution  [2,4,6,12,13,19]. 


The  main  issues  in  the  construction  of  the  Hybrid  method  are  the  smoothness  mea¬ 
surement  of  the  solution  and  the  subdomains  types  switching  algorithm.  In  this  work 
we  employ  the  high  order  multi-resolution  algorithm  by  Ami  Harten  [16]  to  build  a  lo- 
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cal  classification  of  the  solution  into  smooth  and  rough.  Originally  built  to  decrease  the 
work  of  the  fluxes  computations  of  Conservation  Laws,  Harten’s  Algorithm  proposes 
to  use  information  from  coarser  grids  when  the  solution  is  locally  over-represented. 
We  instead  use  the  multi-resolution  information  to  apply  distinct  numerical  method¬ 
ologies  to  the  different  structures  of  the  solution.  The  main  goal  is  to  conjugate  the 
higher  efficiency  of  the  spectral  method  with  the  shock-capturing  capability  of  the 
WENO  method.  The  multi- resolution  analysis  is  used  to  trigger  the  switching  algo¬ 
rithm  to  change  the  subdomains  spatial  discretizations  if  shocks  start  to  develop  at 
a  spectral  subdomain,  or  if  the  solution  becomes  smooth  at  a  WENO  one.  Moving 
discontinuities  are  similarly  treated  by  changing  to  (or  maintaining  as)  WENO  the 
subdomains  on  their  paths  and  switching  to  (or  maintaining  as)  spectral  the  subdo¬ 
mains  that  were  left  behind.  These  changes  are  performed  via  Lagrangian  and  spectral 
interpolations  of  the  local  solutions  to  the  new  discretizations  grids.  Interpolation  is 
also  used  to  patch  the  solutions  at  the  interfaces.  While  a  simple  average  is  sufficient 
for  the  interfaces  where  the  solution  is  smooth,  using  the  same  grid  spacing  at  adjacent 
WENO  subdomains  is  necessary  for  a  conservative  transmission  of  shocks  [28].  Even 
though  we  do  not  have  a  theoretical  proof  of  the  conservation  of  the  Hybrid  scheme, 
we  argue  that  with  the  conjugation  of  two  conservative  schemes  with  a  conservative 
WENO  interface  and  the  high  order  accuracy  of  the  conservative  spectral  scheme,  the 
conservation  error  should  be  spectrally  small.  We  have  numerically  demonstrated  this 
fact  in  [7]  through  a  long  time  time  integration  of  the  inviscid  Burgers  equations  with 
correct  shock  speed  and  achieved  excellent  agreement  with  the  analytical  solution  of 
the  standard  Riemann  shock-tube  problems,  such  as  the  Lax  and  the  Sod  problem  of 
the  Euler  Equations. 


The  paper  is  organized  as  follows:  Section  6  provides  quick  reviews  on  spectral  and 
WENO  methods.  The  Multi- Resolution  analysis  is  discussed  in  details  at  Section  7 
and  the  Hybrid  Method  is  introduced  at  Section  8.  The  Switching  Algorithm  is  pre¬ 
sented  in  Section  8.3  and  numerical  experiments  with  two  dimensional  compressible 
flows  are  finally  presented  at  Section  9.  Concluding  remarks  are  given  in  Section  4. 
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6  Spectral  and  WENO  Methods 


6.1  The  Spectral  Method 


In  spectral  collocation  methods,  the  function  u(x)  is  interpolated  by  a  global  La- 
grangian  interpolation  polynomial  of  degree  N,lj(x),  at  a  given  set  of  collocation 
points  {xj,j  =  0, . . . ,  N}  as 


N 

INu(x )  =  (31) 

j=o 

where  IN  is  the  Interpolating  operator  and  lj{xj)  =  Sij. 


The  Lagrangian  interpolation  polynomial  lj(x)  can  be  constructed  as 


lj(x) 


q(x) 

{x-x^q'ixjY 


N 

q(x)  =  n 

3=0 


and  the  derivative  of  INu  becomes 


(32) 


dINu(x) 

dx 


N 

Y,  u(xi) 


3=0 


dlj  (x) 
dx 


(33) 


In  this  study,  we  will  employ  the  Chebyshev-Gauss-Lobatto  quadrature  points,  namely, 


Xj 


(34) 


which  are  the  roots  of  (1  —  x2)T'N(x )  and  TN(x)  is  the  A^_th  degree  Chebyshev  poly¬ 
nomial  of  the  first  kind.  The  Chebyshev  interpolation  polynomial  is  given  by 


CjN2(x  —  Xj ) 


(35) 
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where  Cj  —  1,  j  —  1, . . . ,  N  ■*=•  1  and  c0  =  cN  =  2. 


The  form  of  the  differentiation  matrix  D  =  dlj^ 1  -  in  (33)  can  be  found  in  [8].  The 
mapping  devised  by  Kosloff  and  Tal-Ezer  [22]  will  also  be  employed  to  enhance  the 
stability  of  the  Chebyshev  collocation  scheme  [5,9]. 


6.2  Filters 


Spectral  methods  are  highly  efficient  and  accurate,  when  the  solution  and  its  deriva¬ 
tives  are  smooth.  In  the  presence  of  discontinuities,  however,  Gibbs  Phenomenon 
generates  oscillations  that  contaminate  the  solution  and  causes  the  loss  of  the  expo¬ 
nential  convergence.  This  is  explained  in  spectral  space  by  the  linear  decay  rate  of 
the  coefficients  an  of  the  global  expansion: 

OO 

U(x)  =  J2anTn(X)-  (36) 

n= 0 

In  such  a  situation,  one  can  modify  the  global  expansion  coefficients  an  to  enhance 
the  convergence  properties  of  the  approximation  via  a  filter  function  <7(77)  [30]  with 
the  following  properties 


<r(v)  =  <r{-v),  o-(±1)  =  0, 

cr(0)  —  1,  a(9)(0)  =  0  q  =  l,...,p-l.  (37) 

If  <7(77)  has  at  least  p  —  1  continuous  derivatives,  <7(77)  is  termed  a  p_th  order  filter.  It 
was  proved  in  [30]  that  the  filtered  expansion  converges  faster  to  the  correct  solution 
than  the  unfiltered  original  one  in  the  case  of  a  discontinuous  function.  Moreover, 
the  convergence  rate  depends  solely  on  the  order  and  the  compactness  of  the  filter 
function  cr (77)  and  the  distance  from  the  discontinuities. 


The  filter  function  used  in  this  study  is  the  Exponential  hlter  given  by 

<7(77)  =  exp(— arf),  (38) 

where  a  =  —  ln(e)  and  e  is  the  machine  zero. 


The  spectral  filtering  can  be  expressed 
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•  in  the  physical  space  as: 


N 

FNu(xi)  =  J2u(xj)$ji,  $ji 
j= 0 


(39) 


•  in  the  transformed  space  as: 


N 

FNu(xi)  =  Y  anc 

n= 0 


T 

-*■  II 


Xj 


2 

Ncn 


N  1 

ujTn(xj)i 

j= 0  ci 


(40) 


where  Fn  is  the  filtering  operator. 


While  filtering  techniques  can  improve  the  overall  convergence  properties  of  the  so¬ 
lution  away  from  discontinuities,  the  solution  near  the  discontinuities  remains  poor. 
Post-processing  of  the  resulting  oscillatory  data  to  recover  spectrally  accurate  non- 
oscillatory  results  can  be  performed  by  various  reconstruction  techniques,  such  as  the 
direct  and  inverse  Gegenbauer  reconstruction  [15,20]  and  Pade  reconstruction  [25]. 
Reconstruction  techniques  are,  in  general,  computationally  costly  and  certain  com¬ 
plications  might  arise.  For  instance,  as  the  degree  of  the  reconstructed  polynomial 
increases,  the  Gegenbauer  transformation  matrices  become  ill-conditioned  clue  to  the 
round-off  error  [20].  Moreover,  the  extension  of  these  reconstruction  techniques  to 
higher  dimensions  is  not  trivial. 


6.3  The  WENO  Method 


WENO  schemes  were  designed  for  the  numerical  solution  of  ffyperbolic  Conserva¬ 
tion  Laws  in  the  form  (for  ease  of  presentation  the  discussion  is  based  on  the  one¬ 
dimensional  formulation) 

ut  +  f(u)a.  =  0.  (41) 

The  Jacobian  Jh  for  (41)  js  needed  to  project  the  system  into  its  characteristic  form 
and,  in  general,  is  not  constant.  To  overcome  the  difficulty  when  computing  the 
flux  at  a  cell  boundary  x _  4 ,  ’’local  freezing”  of  the  matrix  components  is  employed 

using  the  Roe  average  u,  1  [27]  defined  as 

i+2 

/(Ui+i)  -  /(Ui)  =  f'(u  1  )(ui+1  -  Ui).  (42) 
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The  numerical  scheme  for  (41)  can  be  written  in  the  conservative  form  as 


duiit ) 

dt 


1 

Ax 


(43) 


where  is  the  numerical  approximation  to  the  point  value  u(xi,t )  and  Ax  is  the 
uniform  grid  spacing.  The  numerical  flux  based  on  the  {uj,j  =  i  —  r, . . . ,  i  +  s} 


/,,1  =  f(ui-r,...,'Ui+a),  (44) 

satisfies  the  following  conditions: 


•  /  is  Lipschitz  continuous  in  all  arguments; 

•  /  is  consistent  with  the  physical  flux  /,  i.e.  f(u, . . .  ,u)  =  f(u). 

The  solution  to  the  conservative  scheme  (43),  if  converges,  will  converge  to  a  weak 
solution.  This  is  known  as  the  Lax-Wendroff  theorem. 


The  earlier  works  of  van  Leer  [29],  Boris  and  Book  [1]  and  Harten’s  TVD  schemes 
have  led  to  the  introduction  of  Essentially  Non- Oscillatory  (ENO)  schemes  [18],  where 
its  basic  premise  is  that  of  adaptivity  of  the  stencil,  based  on  the  local  smoothness 
of  the  solution,  using  only  a  stencil  free  of  discontinuities  in  the  computation  of  flux 
gradients.  ENO  schemes  have  been  shown  to  generate  sharp  resolutions  of  shocks  as 
well  as  to  maintain  high  order  of  accuracy  in  smooth  regions. 


An  improvement  of  finite  volume  ENO  -  Weighted  Essentially  Non- Oscillatory  (WENO) 
schemes  were  proposed  by  Liu,  Osher,  and  Chan  [23].  WENO  schemes  use  the  same 
idea  as  ENO,  except  that  WENO  uses  a  convex  combination  of  all  available  smooth 
stencils  to  obtain  higher  order  of  accuracy  than  the  original  ENO  at  smooth  parts  of 
the  solution.  Finite  Difference  WENO  schemes  were  later  proposed  by  Jiang  and  Shu 

[19]- 


In  all  numerical  examples  that  follows,  we  use  the  fifth  order  characteristic-wise 
WENO  finite  difference  formulation,  making  use  of  the  Roe  average  for  the  eigen- 
system  defined  above  and  the  global  Lax-Fredrichs  flux  splitting 

/±  =  £(/(«)  ±  au),  (45) 

where  a  =  maxu  max!<j<jvA  |Aj(w)|  and  A i(u),i  =  1, . . .  ,N\  are  the  local  eigenvalues 

of  the  Jacobian  4A  The  numerical  flux  /.  i  is  computed  by  the  WENO  reconstruction 
U  i+2 
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procedure.  The  interested  reader  is  referred  to  [4]  for  details. 


7  Multi-Resolution  Analysis 


The  successful  implementation  of  the  Hybrid  method  depends  on  the  ability  to  obtain 
accurate  information  on  the  smoothness  of  a  function.  In  this  work,  we  employ  the 
Multi- Resolution  (MR)  algorithms  by  Harten  [16,17]  to  detect  the  smooth  and  rough 
parts  of  the  numerical  solution.  The  general  idea  is  to  generate  a  coarser  grid  of 
averages  of  the  point  values  of  a  function  and  measure  the  differences  (MR  coefficients) 
di  between  the  interpolated  values  from  this  sub-grid  and  the  point  values  themselves. 
A  tolerance  parameter  eMR  is  chosen  in  order  to  classify  as  smooth  those  parts  of  the 
function  that  can  be  well  interpolated  by  the  averaged  function  and  as  rough  those 
where  the  differences  di  are  larger  than  the  parameter  eMR.  We  shall  see  that  the  order 
of  interpolation  is  relevant  and  the  ratio  between  di  of  distinct  orders  may  also  be 
taken  as  an  indication  of  smoothness. 


Let  us  start  by  showing  two  examples  where  one  can  notice  the  detection  capabilities 
of  the  Multi-Resolution  analysis  that  will  be  presented  below.  The  left  and  right 
figures  of  Figure  10  show  the  piecewise  analytic  function 


f(x)  =  < 


10  +  x3  —l<x<  —0.5 
x3  —0.5  <  x  <  0 
sin(27nr)  0  <  x  <  1 


(46) 


and  the  density  (p)  of  the  Mach  3  Shock-Entropy  wave  interaction  problem  [19]  as 
computed  by  the  classical  fifth  order  WENO  finite  difference  scheme,  respectively. 


The  test  function  (46)  has  a  jump  discontinuity  at  x  —  —0.5  and  a  discontinuity  at  its 
first  derivative  at  x  —  0.  One  can  see  that  at  each  grid  point  the  differences  dt  decay 
exponentially  to  zero  inside  the  analytical  pieces  of  the  function  when  the  order  of 
interpolation  increases  from  nMR  =  3  to  nMR  =  8.  At  the  discontinuity  x  =  0.5,  the 
measured  differences  di  are  0(1)  and  remain  unchanged  despite  the  increase  of  the 
interpolation  order.  Similar  behavior  is  exhibited  at  the  derivative  discontinuity  at 
x  =  0  with  a  smaller  amplitude. 


Also,  in  the  right  figure  of  Figure  10,  the  density  of  the  Mach  3  Shock-Entropy  wave 
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tion.  (Right)  The  third,  fifth  and  seventh  order  MR  coefficients  di  of  the  density  f(x)  =  p 
of  the  Mach  3  Shock-Entropy  wave  interaction  problem. 

interaction  problem  and  the  corresponding  MR  coefficients  di  are  shown  for  the  third, 
fifth  and  seventh  order  Multi-Resolution  analysis.  The  location  of  the  main  shock  is 
at  x  ~  2.73  and  the  shocklets  behind  the  main  shock  are  well  captured.  The  high 
frequencies  behind  the  main  shock  are  much  better  distinguished  with  the  higher 
orders. 


Averaging  a  function  corresponds  to  filter  the  upper  half  of  the  spectrum.  The  main 
idea  of  Harten’s  smoothness  classification  is  to  measure  how  distant  the  actual  values 
of  the  function  are  from  being  predicted  through  interpolation  of  the  lower  half  of  the 
frequencies  contained  in  the  sub-grid  of  averages.  We  now  describe  a  detailed  con¬ 
struction  of  the  sub-grid  of  averages  and  its  corresponding  interpolating  polynomial, 
finishing  with  a  worked  example. 


Given  an  initial  number  of  grid  points  No  and  grid  spacing  Axo,  consider  the  set  of 
nested  dyadic  grids  { Gk ,  0  <  k  <  L},  defined  as: 


Gk  =  {x^i  =  0,...,Nk},  (47) 

where  xk  =  iAxk,  Ax*,  =  2kAxo,  Nk  =  2~kN0.  For  each  level  k  >  0  we  define  the  set 
of  cell  averages  {fk,  i  =  1 , ,Nk}  at  xk  of  a  function  /(x): 


fl 


1 

A  xk 


f(x)dx, 


(48) 


and  /)°  =  /°.  Let  fki_l  be  the  approximation  to  fki_\  by  the  unique  polynomial  of 
degree  2s  that  interpolates  fk+l,  |Z|  <  s  at  xk+l,  where  r  =  2s  +  1  is  the  order  of 
approximation. 
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The  approximation  differences,  also  called  multiresolution  coefficients,  d/-  =  f^_ \  — 
f-2il i,  at  the  &_th  grid  level  and  grid  point  Xi ,  have  the  property  that  if  f(x )  has  p—  1 
continuous  derivatives  and  a  jump  discontinuity  at  its  p_th  derivative,  then 

dk  _  |  M-[//P)]  for 
[  AxlfP  for 

where  [•]  denotes  the  magnitude  of  the  jump  of  the  function  inside. 


p  <  r 
p  >  r 


(49) 


From  formula  (49)  it  follows  that 

l°^Lrl|  ~  2~P\di\,  where  p  =  min{p, r}.  (50) 


Equation  (50)  shows  that  away  from  discontinuities,  the  MR  coefficients  d\  diminish 
in  size  with  the  refinement  of  the  grid;  close  to  discontinuities,  they  remain  the  same 
size,  independent  of  k.  The  MR  coefficients  df  were  used  in  [17]  in  two  ways.  First, 
finer  grid  data  /)°  were  mapped  to  its  M  level  multiresolution  representation  ff  = 
(d] ,  •  •  •  ,  df,  f\M)  to  form  a  multiscale  version  of  a  particular  scheme,  where  truncation 
of  small  quantities  with  respect  to  a  tolerance  parameter  decreased  the  number  of 
flux  computations.  Secondly,  the  MR  coefficients  df  also  acted  as  a  shock  detection 
mechanism  and  an  adaptive  method  was  designed  where  a  second-order  Lax-Wendroff 
scheme  was  locally  switched  to  a  first-order  accurate  TVD  Roe  scheme,  whenever  d\ 
was  bigger  than  eMR. 


Equation  (49)  also  indicates  that  the  variation  of  the  MR  order,  nMR ,  can  give  ad¬ 
ditional  information  on  the  type  of  the  discontinuity.  Nevertheless,  in  this  work,  we 
will  be  limited  at  using  only  the  first  level  k  —  1  of  the  multiresolution  coefficients 
and  we  shall  drop  the  superscript  1  from  the  d.j  from  here  on  unless  noted  otherwise. 


Hence,  to  find  di,  the  idea  is  to  construct  a  piecewise  polynomial  Pk(x )  of  degree 
k  =  nMR  using  k  + 1  computed  average  values  of  /),  /*,  at  the  equi-spaced  grid  x^  such 
that 


Pk(xi)  =  f(x i)  +  0(  Axk+1), 


(51) 


and 


fi  Pkip^i)' 


(52) 


Given  a  tolerance  level  eMR,  the  smoothness  of  the  function  f(x)  at  Xi  would  then  be 
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checked  against  the  magnitude  of  the  d,,  namely: 


!|rfj|  <  eMR  solution  is  smooth. 

\di\  >  eMR  solution  is  non-smooth. 

The  algorithm  for  computing  the  MR  coefficients  d*  is  given  next. 


(53) 


7.1  Computing  the  MR  Coefficients 


Consider  an  equi-spaced  grid  {xi  =  iAx,  i  =  —  m, . . . ,  0, . . . ,  N, . . . ,  N  +  M }  where 
Aa:  is  the  constant  grid  spacing.  N  can  be  an  odd  or  an  even  number.  Depending  on 
N  and  the  even  or  odd  order  of  the  MR  Analysis  nMR)  the  number  of  ghost  points  m 
and  M  required  are  given  in  table  VIII. 


N 

Umr 

m 

M 

odd 

odd 

n.MR  +  1 

n>MR  +  1 

even 

odd 

UmR  +  1 

UmR 

odd 

even 

UmR 

Umr  +  2 

even 

even 

UmR 

n>MR  +  1 

Table  VIII 

The  number  of  ghost  points  m  and  M  required  for  the  MR  Analysis. 


Given  the  grid  point  values  of  the  function  f(x),  the  average  values  are  computed  as 


fi 


(/A  +  /W  i ) , 


m  N  +  M  —  1 
~2’  7  2 


(54) 


We  construct  a  piecewise  k  =  nMR  degree  polynomial  Pfc(x)  using  the  k  +  1  computed 
average  values  of  the  given  function,  fi  such  that 


Pk(x  i)  =  f(x i)  +  0(Axt+1). 


(55) 


The  polynomial  Pk(xi ),  /  =  \m  and  L  —  l  —  1  or  L  —  l  if  k  is  odd  or  even,  respectively, 
can  be  written  as 

l-\-L 

Pk{xi)  =  J2  arfr.  (56) 

r=i—l 
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However,  since  the  coefficients  a  depend  only  on  xt  and  do  not  depend  on  the  function 
fix),  the  Pk(Xi)  can  be  written  as: 


Pk(xf) 


J2r=-l  Otrfi+r,  Hlod(i,  2)  =  0 

Er=-jA-+l/i+r>  Hlod(i,  2)  =  1 


In  the  case  of  mod(i,  2)  =  1, 


(57) 


/3-r  =  ar ,  r  =  — l , . . . ,  L. 


(58) 


Furthermore,  if  nMR  is  even,  the  coefficients  a  are  symmetric  about  r  =  0,  namely, 
OL—r  =  ar,  r  —  1, . . . ,  L. 


The  desired  coefficients  a  are  computed  by  requiring  Pk{x )  to  be  equal  to  each  of  the 
first  k  +  1  monomials  f(x )  =  1,  x,  x2, . . . ,  xk  and  evaluated  at  any  grid  point  x  =  x*. 
For  simplicity,  we  take  x*  =  0.  The  ft  are  evaluated  for  i  =  — /, . . . ,  L.  This  procedure 
results  in  a  system  of  linear  equations,  A  a  =  b,  where 


f  1 

1  1 

1 

M 

A  = 

-21  +  (-2/  +  1)  .. 

. .  2L+(2L  +  1) 

,  ~oi  = 

01—1+ 1 

,~b  = 

0 

K(-2l)k  +  (-21  +  l)k  .. 

..  (2L)k  +  (2L  +  l)fc  j 

l  ) 

W 

(59) 


and  A  is  a  matrix  of  size  (L  +  l  +  1)  x  {L  +  l  +  1). 


Using  (57),  the  k-th  order  Multi- Resolution  coefficients  d*  at  x,  can  be  computed  as 


di  =  fi-  Pkixi)  i  =  0, . . . ,  N. 


(60) 


One  can  also  evaluate  the  a  by  matching  the  terms  in  the  Taylor  series  expansion  using 
(55)  and  (56)  to  any  desired  order,  however  this  procedure  may  become  cumbersome 
for  high  order  k. 
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Example 


To  illustrate  the  procedure  above,  we  will  construct  two  unique  local  polynomials  with 
k  =  nMR  =  3,  such  that  Pk(x 0)  =  f(x0 )  +  0(Aa;fc+1)  and  Pk(%i)  =  f(x i)  +  0(Aa;fc+1). 


To  construct  the  desired  polynomials  one  needs  to  find  the  unique  coefficients  {a_2,  a_i,  a0,  aq} 
and  p0,  Pi,  /32}  such  that 


a_ 2/_2  +  a_i/_i  +  a0f0  +  «i/i  —  /  (xq)  +  0(  Ax4) 


and 

P-if-i  +  A)/o  +  /3i/i  +  ^2/2  =  /(^l)  +  0(Aa;4). 
The  system  of  equations,  (59),  becomes 


1  1  111^ 

M 

M 

-7-3  15 

(X—l 

0 

,  ~a  = 

,  b  = 

25  5  1  13 

OL  0 

0 

^  —91  —9  1  35  y 

l  ) 

W 

Solving  this  system  yields 


ot- 2  —  — 


64’ 


Ql  —  1  — 


17 

64’ 


«o  — 


55 

64’ 


Qlx  —  — 


64’ 


and  {/3_i  =  cm,  =  aQ,  =  a_i,/32  =  «-2}- 


(61) 

(62) 


(63) 


(64) 


The  tolerance  parameter  eMR  determines  the  dynamic  activation  of  the  spectral  and 
WENO  spatial  discretizations  along  the  various  subdomains  of  the  hybrid  method. 
While  a  too  small  value  of  eMR  activates  the  more  expensive  WENO  method  at  sub- 
domains  where  the  solution  is  smooth,  a  larger  value  activates  the  spectral  method 
at  a  subdomain  with  low  spatial  resolution,  generating  oscillations.  eMR  also  bears  a 
straight  relation  with  the  interpolation  order  nMR.  High  nMR  values  decrease  the  size 
of  eMR  one  needs  to  chose,  since  high  frequencies  are  less  mistaken  by  gradient  jumps. 
The  general  guideline  is  to  start  with  a  value  for  nMR  at  least  equal  to  the  order  of 
the  WENO  method  and  increase  it  according  to  the  complexity  of  the  solution.  For 
instance,  nMR  =  5  is  a  good  choice  for  the  piecewise  smooth  solution  of  the  SOD 
problem,  the  Entropy  problem  would  work  better  with  nMR  =  7.  For  most  of  the 
flows  with  shock  that  were  tested,  the  value  of  eMR  =  1CT3  yielded  a  good  balance 
between  computational  speed  and  accuracy  of  the  numerical  solution. 
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11.  (Le: 

1) 

-‘artitioi 

(Middle)  Typical  spectral  subdomain.  (Right)  Typical  WENO  subdomain. 

8  The  Multi-Domain  Hybrid  Spectral- WENO  method 


We  now  describe  the  implementation  of  the  Multi-Domain  Hybrid  Spectral- WENO 
method  (Hybrid),  detailing  the  structure  of  the  two-dimensional  grid  of  subdomains, 
the  interfaces  treatment  and  the  algorithm  that  switches  the  spatial  discretization  of 
the  subdomains. 


The  main  idea  of  the  Hybrid  method  can  be  formulated  as: 

Avoid  Gibbs  phenomenon  by  keeping  discontinuities  at  WENO  subdomains  and  in¬ 
crease  the  numerical  efficiency  by  treating  the  smooth  parts  of  the  solution  using  a 
spectral  spatial  discretization. 


In  this  study,  the  physical  domain  is  restricted  to  rectangular  shapes  and  will  be 
partitioned  into  a  (Nf  x  Nf)  grid  of  subdomains.  Figure  11  shows  an  example  of  such 
a  domain  partition  along  with  typical  spectral  and  WENO  subdomains.  Note  that 
patching  of  subdomains  occurs  due  to  the  ghost  points  of  the  WENO  discretization. 


We  shall  use  a  vector  k  =  (kx,ky):kx  =  1, . . . ,  Nf,  ky  =  1  to  denote  the 

coordinates  of  the  two  dimensional  subdomain  grid.  For  example,  k  =  (2,  3)  means 
the  subdomain  number  kx  =  2  in  the  x  direction  and  ky  =  3  in  the  y  direction.  Each 
subdomain  is  initialized  either  as  a  spectral  or  WENO  subdomain.  The  Chebyshev- 
Gauss-Lobatto  points  are  used  for  the  spectral  discretization  and  an  uniformly  spaced 
grid  with  ghost  points  is  used  as  the  WENO  grid. 


For  the  sake  of  simplicity,  we  consider  only  square  subdomains,  with  Ns  x  Ns  grid 
points  at  a  pectral  subdomain  and  Nw  x  Nw  equi-spaced  points  at  a  WENO  sub  do- 
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Fig.  12.  (Left)  Typical  WENO  subdomain.  (Middle)  ’’Ghost  Area”.  (Right)  ’’Buffer  Area”. 

main.  The  number  of  ghost  points  is  denoted  by  r  and  NB  is  the  number  of  points 
used  for  the  buffer  zone  (see  below). 


8.1  Description  of  WENO  subdomains 


Each  WENO  subdomain  k  is  composed  of  three  parts:  The  ’’Ghost  Area”,  the  ’’Buffer 
Area”,  and  the  ’’Interior  Area”  as  shown  in  Figure  12. 


•  Ghost  Area: 

The  ’’Ghost  Area”  is  used  for  the  WENO  Reconstruction  and  for  communication 
with  its  neighboring  subdomains  and  is  subdivided  into  eight  ’’Ghost  Zones” 

{G^,  i  —  1, . . . ,  8}  (see  the  middle  figure  of  Figure  12), 


G%  =  If  x  Jf  ,  G*  =  If  x  Jf  ,  =  If  x  Jf 

<  Gf  =  If  x  Jf  ,  ,  =  If  x  Jf  , 

G*2  =  If  x  Jf  ,  G*3  =  If  X  Jf  ,  G\  =  If  X  Jf 


(65) 


where  the  stencil  ranges  are 


If  —  { T  — r,  •  •  •  1  ^-l})  If  —  {To>  ■  ■  ■  )  XN},  If  —  {Xjv+1,  •  •  •  ,  iCjv+r}) 

Iq  {V—rt  •  •  •  1  V—  1 } )  {jo  1  •  •  ■  1  Vn  } )  J2  {Vn+1  )  •  •  •  )  Vn+t}  ) 


with  N  =  Nw. 

•  Buffer  Area: 

The  Buffer  Area  is  used  for  the  treatment  of  moving  discontinuities.  If  a  high 
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gradient  or  shock  is  detected  inside  a  buffer  zone,  the  closest  neighboring  spectral 
subdomain(s)  is  switched  to  WENO  subdomain(s). 

For  two-dimensional  problems,  each  WENO  subdomain  has  eight  ’'Buffer  Zones” 
{Bf,  i  =  A, ... ,  H}  (see  the  right  figure  of  Figure  12), 


Bl 

<  B\ 

Bl 


Ib0  x  ,  Bk  =  I\  x  J\  ,  B k  =  I\  x  J\ 

Ib0  x  Jb  ,  ,B*  =  Ib  xJb\, 

Ib  x  Jb  ,  Bk  =  If  x  J0b  ,  5k  =  jb  x  jb 


(66) 


where  the  stencil  ranges  are  given  by 


Ib 

1o 

Tb 

J0 


{^o,  ■  ■  ■  i  XM}, 
{Hoi  •  •  •  5  Um },  J\ 


Af+1) 


{^M+l)  •  •  •  )  )  1 2  —  {*£jv- 

•  •  •  )  y N—M  }  1  J2  =  \IJn-M+1i 


with  =  Nw  and  M  =  iVs.  As  greater  is  the  value  of  NB,  earlier  is  the  detection 
of  shocks  and  gradients,  however,  at  the  cost  of  early  switching  of  the  neighboring 
subdomains  to  WENO  and  greater  chance  of  unnecessary  costly  computations. 
Satisfactory  results  have  been  obtained  with  the  default  M  =  r  in  this  study.  It 
should  be  noted  that  these  buffer  zone  grid  points  are  part  of  the  interior  grid 
points  and  should  not  be  confused  with  WENO  ghost  points. 

•  Interior  Area: 

Finally,  the  ” Interior  Area”  is  all  the  WENO  grid  points  Wk  excluding  the  ’’Ghost 
Area”  and  ” Buffer  Area”,  Jk  =  hhk/(Ui  Gk  ©  U*  Bk).  If  the  high  gradient  stays 
inside  this  area,  no  action  is  required  by  the  neighbors. 


8.2  Interface  Conditions 


The  following  configurations  are  representative  of  any  two-dimensional  domain  par¬ 
titions: 


•  Spectral-Spectral-Spectral-Spectral  (Figure  13) 

•  Spectral-Spectral-Spectral- WENO  (Figure  15) 

•  Spectral-Spectral- WENO- WENO  (Figure  15) 

•  Spectral- WENO-WENO-WENO  (Figure  15) 

•  WENO-WENO-WENO-WENO  (Figure  14) 
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The  functional  values  in  the  spectral  and  WENO  subdomains  are  denoted  by  s\:j  and 
wfj  respectively  at  ( x^yj )  in  the  respective  subdomain  k.  Centered  at  the  subdomain 
k,  we  define  the  superscripts  k0,  ki,  k2,  k3,  k4,  k5,  k6,  k7,  k8  as 


k§  =  k  +  ( — 1,  +1)  ,  k7  =  k  +  (+0,  +1)  ,  kg  =  k  +  (+1,  +1) 

ki  =  k  +  ( — 1,  +0)  ,  ko  =  k  ,  k,5  =  k  +  (+1,  +0)  j 

^  k2  =  k  +  (— 1,  —  1)  ,k3  =  k+  (+0,  —1)  ,  k4  =  k  +  (+1,  —  1)  ^ 


(67) 


and  they  denote  the  Center  (k0),  Left  (k4),  Bottom-Left  (k2),  Bottom  (k3),  Bottom- 
Right  (k4),  Right  (k5),  Top-Right  (k6),  Top  (k7),  and  Top- Left (  k8)  subdomains  of 
subdomain  k  =  k0  respectively.  For  example  sff  is  to  be  interpreted  as  a  value  at 
the  point  ( x^yj )  in  the  spectral  Left  subdomain  k4  =  k  +  (— 1,+0).  The  subdo¬ 
main  interfaces  consist  of  corners  and  shared  sides  among  the  spectral  and  WENO 
subdomains. 


Spectral- Spectral  Interface 


1 

1 

— - 

s1 

C 

2 

— 

4 

s2 

1 

S 

s3 

s4 

Fig.  13.  (Left)  Spectral-Spectral-Spectral-Spectral  subdomains  configuration.  (Middle)  Cor¬ 
ner  point  of  all  four  spectral  subdomains.  (Right)  Shared  side  of  two  spectral  subdomains. 


Corners  only  appear  at  the  connection  of  two  or  more  spectral  subdomains,  since  the 
grid  points  of  the  WENO  subdomains  never  coincide  with  the  spectral  collocation 
points  at  the  interface.  Corner  values  are  assigned  with  the  average  of  the  values 
of  all  connecting  subdomains.  For  example,  in  figure  13,  assuming  that  S 1  is  the 
reference  subdomain  k0  and  N  =  Ns, 

e^o  —  gk"3  —  f  e^o  -|-  S^3  _L_  -L  e^3  )  1681 

S0N  '’oo  SN 0  SNN  ^  \°0N  ‘  '’OO  “  & N0  W  O  N NJ  ■  ) 
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Along  the  interface  between  two  spectral  subdomains,  the  values  at  the  shared  side 
are  computed  similarly  by  assigning  the  average  of  the  two  values  of  the  connected 
subdomains  at  all  collocation  points.  Consider  the  interface  between  the  two  spectral 
subdomains  S 1  —  S 2  (the  rightmost  figure  of  Figure  13), 


ko  _  „ks  _ 

0  j  SNj 


j  =  0,...,N. 


(69) 


The  other  spectral  subdomains’  corners  and  sides  are  treated  similarly. 

Note  8.1.  In  case  of  adjacent  subdomains  with  different  numbers  of  collocation  points, 
global  interpolation  can  be  used  before  averaging. 

Note  8.2.  An  exact  Riemannn  solver  and  penalty  interface  conditions  had  been  tried 
but  no  discernible  differences  were  observed  from  the  simple  average. 


WENO-  WENO  Interface 


To  maintain  conservation,  adjacent  WENO  subdomains  are  required  to  have  the  same 
spacing  Ax  (see  [28]).  This  implies  that  ghost  points  of  a  WENO  subdomain  match 
the  interior  points  of  the  neighboring  WENO  subdomain  (see  Figure  14).  Thus,  simple 
copying  of  the  solution  values  of  the  neighboring  interior  points  to  the  ’’Ghost  Area” 
is  sufficient  for  the  communication  of  WENO  subdomains. 


Fig.  14.  (Left)  WENO- WENO- WENO- WENO  subdomain  configuration.  (Middle)  fW-W2 
Interface.  (Right)  W1  —  W 4  Interface. 


For  instance  the  middle  figure  of  Figure  14  shows  two  ’’Ghost  Zones”:  The  zone  in  the 
solid  rectangle  corresponds  to  the  ’’Ghost  Zone”  of  W 1 ;  and  the  zone  in  the  dash- 
dotted  rectangle  corresponds  to  the  ’’Ghost  Zone”  G\  of  W2  (see  65).  In  this  case,  by 
denoting  N  =  Nw  and  using  W1  as  a  reference  subdomain  (subdomain  k  =  k0) 


w 


ko 

(JV+()j 


=  W^5 
wij  i 


w 


ks 

-w 


=  w 


ko 

(N+l-l)j  1 


j  =  -r, . . . ,  N  +  r, 


Z  =  l,...,r,  (70) 
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where  r  is  the  number  of  ghost  points. 

Note  8.3.  In  general  the  number  of  ’’ghost”  points  in  each  subdomain  does  not  have 
to  be  the  same  as  long  as  the  same  grid  spacing  Ax  is  used  at  all  the  neighboring 
subdomains. 


Spectral-  WENO  Interface 


Spectral- WENO  interface  can  be  found  in  the  configurations  shown  in  Figure  15. 


Fig.  15.  Subdomain  Configurations:  (Left)  Spectral-Spectral-Spectral- WENO.  (Middle) 
Spectral-Spectral- WENO- WENO.  (Right)  Spectral- WENO- WENO-WENO. 


The  ’’Ghost  Zones”  G\,  G y  and  G§  of  the  WENO  subdomain  lie  in  the  three  neighbor¬ 
ing  spectral  subdomains  as  shown  in  Figure  15.  The  ’’Ghost  Area”  points  of  WENO 
subdomains  are  computed  via  spectral  interpolations.  To  be  more  specihc,  referring 
to  the  middle  figure  of  Figure  16,  the  solution  in  the  area  of  the  intersection  of  S3 
(subdomain  k)  and  W4  (subdomain  ks)  corresponding  to  the  ’’Ghost  Zone”  in 
(65)  (see  also  Figure  12)  is  read  using  the  interpolating  polynomial  representing  the 
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solution  in  S'3,  that  is,  at  the  subdomain  k, 


NS  NS 

w(xn,  ym)  =  siMxn)lj(ym),  m  =  0, . . . ,  Nw,  n  =  -r, . . . ,  -1,  (xn,  ym)  G  k5, 

j= 0  i=0 

(71) 

where  sVJ  are  the  functional  values  in  S'3  (subdomain  k),  U(x)  and  lj(y)  are  Lagrangian 
interpolation  polynomials  of  degree  Ns.  The  ghost  points  for  the  other  two  configu¬ 
rations  can  be  obtained  in  a  similar  way. 


The  interface  points  of  the  spectral  subdomain  are  computed  via  a  two  dimensional 
interpolation  polynomial  of  degree  r  =  nw.  The  choice  of  stencils  {(xi,yj),i  =  —  i0  — 
r / 2, . . . ,  *0  +  r/2,  j  =  j0  —  r / 2, ... . ,  j0  +  rf  2}  for  the  interpolation  polynomial  should 
be  as  symmetric  about  a  given  point  (xt01  yJ0)  as  possible. 


8.3  The  Switching  Algorithm 


In  this  section  we  describe  the  algorithm  used  to  switch  the  subdomains  spatial  dis¬ 
cretizations  as  indicated  by  the  Multi-Resolution  Analysis  of  Section  7.  The  following 
three  conditions  are  the  main  rules  to  be  followed: 


(1)  If  a  subdomain  contains  high  gradients,  then  switch  its  spatial  discretization  to 
(or  keep  it  with)  WENO; 

(2)  If  high  gradients  are  present  in  the  ’'Buffer  Areas”  of  connected  neighboring 
subdomains,  then  switch  the  current  subdomain  to  (or  keep  it  as)  a  WENO 
subdomain; 

(3)  In  any  other  case,  switch  the  subdomain  to  (or  keep  it  as)  a  spectral  subdomain; 


The  Erst  condition  above  avoids  the  Gibbs  phenomenon,  keeping  the  discontinuities 
inside  WENO  subdomains.  The  second  condition  ensures  the  switch  to  WENO  sub- 
domain  in  order  to  allow  only  WENO-to-WENO  transmission  of  discontinuities.  The 
third  condition  improves  the  numerical  efficiency,  since  it  ensures  that  smooth  parts 
of  the  solution  will  always  be  contained  in  spectral  subdomains. 


Multi-Resolution  analysis  of  the  solution  is  performed  at  every  sub-stage  of  the  third 
order  TVD  Runge-Kutta  scheme  used  for  the  temporal  evolution.  At  each  subdomain 
k,  we  define  the  smoothness  flag  variable,  Flag^,  at  each  grid  point  (ay,  y3)  including 
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the  ghost  points,  as  , 


Flag!" 


i.  Ml5l>£"«  for  (x<’ Vi)  £  ,k ® Bik  (  =  A . H 

0,  otherwise 


(72) 


where  dfj  are  the  MR  coefficients.  Since  the  Multi-Resolution  Analysis  requires  uni¬ 
formly  spaced  grids,  the  spectral  grids  are  first  interpolated  to  uniformly  spaced  grids 
before  obtaining  the  MR  coefficients  dfj.  The  necessary  ghost  points  are  acquired  from 
the  neighboring  subdomains.  At  the  boundary  subdomains,  the  value  of  the  boundary 
ghost  points  are  extrapolated  linearly  from  the  interior  data. 


The  algorithm  proceeds  by  checking  for  each  subdomain  k  and  at  the  buffer  zones  of 
the  neighboring  subdomains  if  Flag^  is  equal  to  one.  If  so,  it  switches  subdomain  k  to 
a  WENO  discretization.  On  the  other  hand,  if  Flag^.  is  identically  zero,  then  a  spectral 
discretization  is  implemented,  or  kept.  These  switches  require  the  use  of  interpolation 
from  a  Chebyshev  grid  of  points  to  a  uniformly  spaced  one  and  vice-versa: 


•  To  switch  from  the  spectral  subdomain  to  the  WENO  subdomain,  the  data  are 
interpolated  onto  the  uniformly  spaced  grid  via  the  spectral  interpolation  formula. 

•  To  switch  from  the  WENO  subdomain  to  the  spectral  subdomain,  the  data  are 
interpolated  onto  the  Chebyshev  Gauss-Lobatto  points  via  the  Lagrangian  inter¬ 
polation  polynomial  of  the  same  order  as  the  WENO  method. 


Back  and  forth  switching  between  WENO  and  spectral  discretizations  may  occur  too 
frequently  for  the  same  domain  when  the  eMR  is  marginally  set.  The  g?F  coefficients 
might  oscillate  around  the  parameter  eMR  in  time  due  to  some  numerical  factors 
such  as  dissipation,  dispersion  and  nonlinear  effects,  or  any  combination  of  such.  This 
pattern  of  switching  can  repeat  itself  for  a  while  until  the  solution  settles  down  with  a 
clear  definition  of  the  dfj,  which  is  either  greater  than  or  smaller  than  the  MR  tolerance 
eMR.  In  order  to  alleviate  such  occurrences,  one  must  devise  a  procedure  preventing 
the  switch  from  WENO  to  spectral  if  it  had  already  occurred  recently.  However, 
the  procedure  must  never  prevent  a  spectral  to  WENO  switch,  for  oscillations  and 
instability  might  occur. 
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9  Numerical  Results 


In  this  section,  we  apply  the  hybrid  method  to  two  well-known  problems  in  Conserva¬ 
tion  Laws:  The  Shock- Vortex  Interaction  and  the  Richtmyer- Meshkov  Instability.  The 
governing  equations  are  the  two-dimensional  Euler  equations  in  Cartesian  coordinates 
given  by: 

Qt  +  Fx  +  Gy  =  0,  (73) 

where 


Q  =  (p,pu,pv,E)T , 

F  =  (pn,  pu 2  +  P,  puv,  (. E  +  P)uj  ,  (74) 

G  =  (pv,  puv ,  pv2  +  P,(E  +  P)v)  , 
and  the  Equation  of  state 

P=(1-l)(E  +  \p(u2  +  v2)),  7  =  1-4.  (75) 


We  will  restrict  the  following  study  to  rectangular  domains,  which  will  be  partitioned 
into  an  equal  number  of  subdomains  in  both  x  and  y  directions.  The  number  of 
Chebyshev  collocation  points  for  all  spectral  subdomains  will  be  the  same  in  both  x 
and  y  directions,  as  well  as  the  number  of  uniformly  spaced  grid  points  for  the  WENO 
subdomains.  The  order  of  the  Multi- Resolution  Analysis  is  the  same  as  the  WENO 
scheme,  nMR  =  nw.  A  14_th  order  Exponential  filter  and  Kosloff-Tal-Ezer  mapping  are 
employed  in  all  spectral  subdomains.  Free  stream  boundary  conditions  are  imposed 
in  the  inflow  and  outflow  in  the  x  direction  and  periodical  boundary  condition  is 
imposed  in  the  y  direction.  To  evolve  the  ODE  from  the  semi-discretized  PDE  in 
time,  the  third  order  Total  Variation  Diminishing  Runge-Kutta  scheme  (RK-TVD) 
will  be  used  [4]: 


U 1  =  Un  +  A  tL(U)n 

V2  =  i(3Pn  +  U1  +AtL(U1))  ,  (76) 

Un+1  =  i  (Un  +  2 U2  +  2A tL(U2)) 

where  L  is  the  spatial  operator.  The  CFL  numbers  for  the  spectral  and  WENO 
subdomains  are  set  to  be  3  and  0.4,  respectively.  All  numerical  experiments  were  run 
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on  a  667  MHz  Compaq  Alpha  machine  with  1GB  memory  and  with  an  Alpha  internal 
floating  point  processor. 


9. 1  Shock  -  Vortex  Interaction 


The  tangential  velocity  profile  of  the  counter-clockwise  rotating  vortex  [21]  centered 
at  (xc,  yc)  and  strength  T  is  given  in  polar  coordinates  by: 


U(r) 


Tr(ro2  —  rf2)  0  <  r  <  ro  <  rq 
<  T r(r-2  —  rf2)  r0  <  r  <  rq 
0  r  >  rq 


where  r0  =  0.2  and  rq  =  1.0. 


(77) 


Due  to  the  strong  nonlinearity  of  the  shock-vortex  interaction  its  physics  are  not 
well  understood.  A  better  understanding  will  have  many  potential  applications,  for 
instance,  subsonic  and  supersonic  jet  nozzle  design  and  blade  noise  generation.  The 
Hybrid  method  will  be  tested  on  this  problem  with  a  vortex  of  amplitude  T  =  0.25 
and  shock  Mach  numbers  Ms  =  1.25,  3,  6. 


Mach  1.25 


In  this  first  example,  the  physical  domain  (0  <  x  <  3.9,  —2  <  y  <  2)  is  partitioned 
into  a  13  x  10  grid  of  subdomains.  Spectral  subdomains  use  a  32  x  32  grid  of  Chebyshev 
points  and  WENO  grids  are  50  x  50.  MR  analysis  is  performed  with  the  MR  tolerance 
set  to  eMR  =  5  x  10~2. 

When  an  initially  planar  shock  wave  hits  the  vortex,  it  deforms  as  compression 
and  rarefaction  regions  are  created  behind  the  shock,  as  shown  in  Figure  17  for 
t  =  0.732, 1.08, 1.2.  As  the  interaction  proceeds  over  time,  strong  bifurcation  and 
deformation  of  the  shock  are  observed.  The  shock  emanating  from  the  compression 
region  has  one  part  moving  upward  and  another  moving  downward.  The  strength  of 
the  lower  upward  moving  part  is  greater,  due  to  the  direction  of  rotation  of  the  vor¬ 
tex.  MR  Analysis  performed  well  in  capturing  the  shock  and  the  high  gradient  regions 
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(a)  t  =  0.73 


(b)  t  =  1.08 


(c)  t  =  1.20 


Fig.  17.  Density  p  of  the  Shock- Vortex  interaction  with  Mach  number  Ms  =  1.25  and 
T  =  0.25  at  (a)  t  =  0.73,  (b)  t  =  1.08  and  (c)  t  =  1.2,  as  computed  by  the  Hybrid  method. 


immediately  behind  the  shock,  as  indicated  by  the  WENO  subdomains  enclosed  with 
black  bounding  boxes.  The  remaining  subdomains  are  accurately  dealt  with  by  the 
spectral  methods.  The  number  of  WENO  subdomains  is  far  fewer  than  spectral  ones 
resulting  in  a  more  efficient  algorithm  than  the  classical  WENO  scheme. 


Mach  3 


We  now  increase  the  Mach  number  of  the  shock  from  Ms  =  1.25  to  Ms  =  3  and 
compare  the  solution  of  the  Hybrid  scheme  with  a  highly  resolved  one  computed  with 
the  classical  fifth  order  WENO  scheme  using  1200  x  1200  points.  The  physical  domain 
(0  <  x  <  3.0,  —2  <  y  <  2)  is  partitioned  into  a  10  x  10  grid  of  subdomains  and  all 
other  parameters  are  as  in  Example  1. 


As  in  the  previous  example,  an  acoustic  wavefront  is  generated  and  a  number  of  fine 
scale  structures  are  formed  behind  the  main  shock.  This  example  indicates  that  the 
Hybrid  method  captures  the  fine  features  of  the  solution  even  in  the  case  of  strong 
shocks. 
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(a)  t  =  0.35,  Hybrid  Scheme  (b)  t  =  0.60,  Hybrid  Scheme 


(c)  t  =  0.35,  WEN05  Scheme  (d)  t  =  0.60,  WEN05  Scheme 

Fig.  18.  Density  p  of  the  Shock- Vortex  interaction  with  Mach  number  Ms  =  3  and  T  =  0.25 
at  time  (a)  t  =  0.35  and  (b)  t  =  0.6  as  computed  by  the  Hybrid  scheme  and  (c)  t  =  0.35 
and  (d)  t  =  0.6  as  computed  by  the  classical  fifth  order  WENO  finite  difference  scheme 
(WEN05)  with  1200  x  1200  grid  points. 


Mach  6 


In  this  last  example,  we  further  increase  the  shock  Mach  number  to  Ms  =  6  and  use 
a  19  x  20  grid  of  subdomains  to  partition  the  physical  domain  (0  <  x  <  4.2,  —3  < 
y  <  3).  The  spectral  grid  is  16  x  16.  This  example  shows  that  the  Hybrid  method 
can  be  applied  to  higher  Mach  number  flows  as  well,  still  resulting  on  a  reliable  and 
efficient  shock-capturing  method. 
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fa)  t  =  0.3 


(b)  t  =  0.4 


Fig.  19.  Density  p  of  the  Shock- Vortex  interaction  with  Mach  number  Ms  =  6  and  T  =  0.25 
at  time  (a)  t  =  0.3  and  (b)  t  =  0.4  as  computed  by  the  Hybrid  method. 

9.2  Two-dimensional  Richtmyer- Meshkov  Instability 


Richtmyer  in  1960  [26]  theoretically  predicted  the  occurrence  of  instability  on  a  per¬ 
turbed  material  interface  under  the  impulsive  acceleration  of  an  incident  shock  wave. 
In  1970  Meshkov  [24]  experimentally  confirmed  these  predictions.  A  variety  of  mo¬ 
tions  can  be  generated  following  the  interaction  of  a  shock  wave  with  an  interface 
separating  two  materials.  Any  small  perturbation  present  on  the  interface  will  be 
amplified  after  such  a  contact.  This  class  of  problems  is  referred  in  the  literature  as 
the  ” Richtmyer-Meshkov  Instability  (RMI)”.  As  the  interface  between  two  materi¬ 
als  becomes  increasingly  distorted  other  instabilities  such  as  the  Kelvin- Helmholtz 
Instabilities  develop  and  a  region  of  turbulence  mixing  ultimately  results.  The  RMI 
arises  in  many  applications  as,  for  instance,  the  Inertial  Confinement  Fusion  (ICF) 
process.  A  recent  model  under  extensive  study  consists  of  a  set  of  laser  beams  directed 
into  a  chamber  containing  a  spherical  fusion  fuel  target.  The  expected  result  should 
be  compression,  ignition  and  a  subsequent  energy  surplus.  However,  since  no  perfect 
capsule  exists,  irregularities  on  the  surface  excite  the  undesirable  RMI,  reducing  the 
effective  uniform  compression  pressure  onto  the  capsule. 
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1) 

2) 

3) 

4) 

5) 

6) 

7) 

8) 

9) 


Rankine-Hugoniot  condition  for  shocks 
Pre-Shock  Temperature  T  =  296  K 
Pre-Shock  Pressure  P  =  0.5  atm 
Xenon  density  pxe  =  2.90  x  10~3 
Argon  density  pAr  =  0.89  x  10~3  ^3- 
Specific  heat  ratio  7  =  | 

Atwood  number  At  =  0.54 
Mach  number  M  =  4.46 
Wave  Length  A  =  3.6cm 


10)  Amplitude  a  =  1.0cm 


Fig.  20.  Initial  Condition  for  the 
tinn 


Richtmyer-Meshkov  Instability  simula- 
Regions  of  Interest  : 

1)  Reflected  shock  generated  by  the  shock 

refraction; 

2)  The  penetration  of  the  heavy  (Xe)  to  light 

(Ar)  gas  forms  the  Spike; 

3)  Triple  point  on  the  transmitted  shock; 

4)  A  small  jet  and  its  vortical  structure; 

5)  The  penetration  of  the  light  (Ar)  to  heavy 

(Xe)  gas  forms  the  Bubble; 

6)  Vortical  rollups  of  the  gaseous  interface. 


Fig.  21.  Typical  regions  of  interest  for  the  simulation  of  the  RMI  at  time  t  =  50  ps.  Only 
the  lower  half  of  the  interface  is  shown. 


Presently,  a  rectangular  domain  with  a  shock  Mach  number  Ms  interacting  with 
a  single  mode  sinusoidal  perturbation  along  a  Xenon  (Xe)  and  Argon  (Ar)  gases 
interface  is  simulated  using  the  Hybrid  method.  The  initial  condition  is  given  in 
Figure  20  and  a  diffusive  interface  is  modeled  with  an  exponential  function,  i.e. 


S(x,y)  =  < 


1  d  <  0 

exp(— ajdl^)  0  <  d  <  1, 

0  d  >  1 


where 


d 


(xs  +  a  cos(27tj//A)  +  5)  —  x 
26 


(78) 


(79) 
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5  =  0.2cm  >  0  is  the  interface  thickness,  (3  =  8  is  the  interface  order,  xs  =  0.5cm 
is  the  location  of  the  interface  and  a  =  —  lne,  where  e  is  the  machine  zero.  The 
conservative  or  primitive  variables  are  scaled  according  to  the  S(x,y)  between  the 
Xenon  and  Argon  gases. 


As  the  shock  wave  collides  with  the  interface  separating  the  two  gases,  the  sine  wave 
perturbation  is  accelerated,  compressed  and  amplified  following  the  refraction  of  the 
shock.  The  heavier  Xenon  gas  (Xe)  will  penetrate  into  the  lighter  Argon  gas  (Ar) 
forming  finger-like  structures  -  bubbles  and  spikes.  A  bubble  (spike)  is  a  portion  of 
the  light  (heavy)  gas  penetrating  into  the  heavy  (light)  gas.  Some  of  these  interesting 
fluid  structures  such  as  the  bubbles,  the  spikes,  the  interfacial  mixing  region  and  the 
vortical  rollup  can  be  observed  at  the  earlier  time  as  depicted  in  Figure  21. 


The  basic  mechanism  of  these  instabilities  is  the  baroclinic  generation  of  vorticity  u 
induced  from  the  misalignment  of  the  pressure  gradient  Vp  of  the  shock  and  the  local 
density  gradient  Vp  across  the  interface: 


dio 

d t 


~  Vp  x  Vp, 


uj  =  V  x  u, 


(80) 


where  u  is  the  velocity. 


Mach  4-46 


In  this  example,  the  physical  domain  (0  <  x  <  5.1,0  <  y  <  A)  is  partitioned  into  a 
17  x  12  grid  of  subdomains.  The  spectral  grids  are  32  x  32  and  the  WENO  grids  are 
50  x  50.  The  MR  tolerance  is  now  lowered  to  eMR  =  5  x  10 ~5.  Once  again,  as  seen 
in  Figure  22,  the  WENO  method  is  activated  only  along  the  material  interface  and 
where  high  gradients  appear. 


Mach  8 


The  physical  domain  (0  <  x  <  24.6,  0  <  y  <  A)  is  partitioned  into  a  82  x  12  grid  of 
subdomains  in  order  to  apply  the  Hybrid  method  for  a  longer  time  integration  up  to 
t  =  lOOps.  The  shock  Mach  number  is  increased  from  Ms  =  4.46  to  Ms  =  8  and  the 
spectral  grid  is  set  to  24  x  24. 
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>•  2 


(c)  t  =  37.5  ps 


(d)  t  =  50.0  ps 


Fig.  22.  Contour  plot  of  the  density  p  with  Ms  =  4.46,  A  =  3.6cm,  a  =  1cm  at  time  (a) 
t  =  12.5  p,s,  (b)  t  =  25.0  ps,  (c)  t  =  37.5  ps  and  (d)  t  =  50.0  ps  of  the  Richtmyer-Meshkov 
Instability  as  computed  by  the  Hybrid  scheme. 


It  can  be  observed  in  Figure  23  that  the  Hybrid  method  successfully  tracks  shocks  and 
high  gradients  with  WENO  discretization  (black  bounding  boxes)  while  the  smooth 
parts  of  the  solution  are  well  represented  using  spectral  subdomains. 
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(a)  t  =  100.0  ps 


(b)  t  =  100.0  ps  with  zoom 


Fig.  23.  Contour  plot  of  the  density  p  with  Ms  =  8,  A  =  3.6cm,  a  =  lcm  at  time  (a) 
t  =  18.5  ps,  (b)  t  =  52.0  ps,  (c)  t  =  100.0  ps  and  (d)  t  =  100.0  ps  with  zoom  of  the 
Richtmyer-Meshkov  Instability  as  computed  by  the  Hybrid  scheme. 


9.3  CPU  Timing 


The  Hybrid  method  has  the  potential  advantage  of  being  faster  than  the  classical 
WENO  method  due  to  the  higher  numerical  efficiency  of  spectral  methods  at  smooth 
parts  of  the  solution.  Moreover,  spectral  discretizations  avoid  the  expensive  character¬ 
istic  decompositions  and  projections  of  the  WENO  method.  In  this  section,  we  provide 
some  CPU  timing  results  for  the  Mach  3  Shock- Vortex  Interaction  (see  section  9.1) 
when  using  the  Hybrid  and  the  classical  WENO  methods  with  equal  resolution. 


57 


The  physical  domain  is  partitioned  into  a  set  of  subdomains  of  sizes  ((10  x  2J)  x  (10  x 
2 J)),j  =  0, 1,2,3.  Spectral  subdomains  use  a  12  x  12  grid  of  Chebyshev  points  and 
WENO  ones  use  20  x  20  grids  of  uniformly  spaced  points.  The  classical  fifth-order 
WENO  method,  here  denoted  WEN05,  uses  the  corresponding  number  of  grid  points 
as  if  all  the  subdomains  in  the  Hybrid  method  were  WENO  subdomains.  The  MR 
tolerance  used  is  eMR  =  5  x  10-2. 


Number  of  subdomains 

Grid  size 

Hybrid  S12W20 

WENO  5 

Speedup 

10x10 

200x200 

265 

282 

1.06 

20x20 

400x400 

2009 

2762 

1.37 

40x40 

800x800 

14410 

26090 

1.81 

80x80 

1600x1600 

112900 

253996 

2.24 

Table  IX 

CPU  timing  in  seconds  and  speedup  factor  for  the  Shock- Vortex  problem  at  time  t  =  0.6 
as  computed  by  the  Hybrid  (with  constant  eMR  =  5  x  10-2)  and  the  WEN05  methods. 


Table  IX  shows  that  a  significant  speed  up  is  achieved  when  using  the  Hybrid  method 
over  the  classical  WEN05  for  increasing  resolution.  Figure  24  shows  the  history  of 
the  coverage  of  the  WENO  subdomains  as  a  percentage  of  the  total  number  of  the 
subdomains.  This  percentage  varies  between  10%  —  20%  for  the  10  x  10  subdomain 
partition  and  gets  proportionally  smaller  by  a  factor  of  2  when  the  number  of  subdo¬ 
main  partition  increases  by  a  factor  of  2  in  each  direction. 
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Fig.  24.  The  time  history  of  the  coverage  of  the  WENO  subdomains  as  a  percentage 
of  total  number  of  subdomains  for  the  four  subdomain  partitions:  (from  left  to  right) 
(10  x  10),  (20  x  20),  (40  x  40),  (80  x  80). 


Reducing  the  MR  Tolerance  level  of  last  experiment  to  eMR  =  5  x  10~3,  while  keeping 
all  the  other  parameters  fixed,  leads  to  a  slight  increase  on  the  CPU  time  usage  for 
the  Hybrid  method  as  compared  with  the  previous  case,  as  shown  in  Table  X.  With 
a  lower  MR  tolerance,  more  subdomains  are  classified  as  containing  high  gradients. 


Number  of  subdomains 

Grid  size 

Hybrid  S12W20 

WENO 

Speedup 

10x10 

200x200 

332 

282 

none 

20x20 

400x400 

2239 

2762 

1.23 

40x40 

800x800 

15580 

26090 

1.67 

Table  X 

CPU  timing  in  seconds  and  speedup  factor  for  the  Shock- Vortex  problem  at  time  t  =  0.6 
as  computed  by  the  Hybrid  (constant  eMR  =  5  x  10~3)  and  WEN05  methods. 
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10  High  Order  Simulation  of  the  Richtmyer-Meshkov  Instability 


We  also,  in  collaboration  with  Dr.  Oleg  Schilling  at  the  Lawrence  Livermore  National 
Laboratory  (LLNL)  pursuing  the  investigation  of  the  high  order  WENO  methods 
for  the  three  dimensional  Richtmyer-Meshkov  Instability.  We  have  recently  extended 
the  Euler  Equation  to  the  full  Navier-Stokes  equation  with  multiple  species  in  two 
and  three  dimensions.  We  are  also  in  the  process  of  implementing  a  full  turbulence 
analysis  routines  for  real  time  data  extraction  and  analysis. 


A  systematical  and  self-consistent  study  was  conducted  to  validate  the  WENO  meth¬ 
ods  for  the  RMI  problems  with  reshock  against  the  Mach  1.21  shock  tube  experiments 
of  Collins  and  Jacobs.  The  qualitative  comparison  above  shows  that  it  is  possible 
to  achieve  very  good  agreement  between  a  two-dimensional,  high-resolution  shock- 
capturing  simulation  with  high-order  flux  reconstruction  and  experimental  density 
PLIF  images.  The  simulations  at  5  and  6  ms  demonstrates  that  higher-order  recon¬ 
struction  better  captures  secondary  instabilities  on  the  interface  and  within  the  vortex 
cores.  The  reshock  takes  place  at  about  t  =  5.75  ms.  In  addition,  the  roll-ups  in  the 
simulation  appear  tighter  and  sharper,  and  more  fine-scale  structures  are  present. 
The  fields  from  the  numerical  simulations  also  supplement  the  experimental  images 
by  displaying  the  shock  focusing  observed  during  the  reshock  process.  This  results  in 
the  formation  and  persistence  of  large-scale  structures  in  the  simulations,  consistent 
with  the  inverse  cascade  of  kinetic  energy  from  small  scales  to  larger  scales  observed 
in  experiments  and  simulations  of  two-dimensional  turbulence. 


As  the  results  of  this  study,  we  have  published  three  papers  in  the  Physics  of  Fluids, 
Journal  of  Computational  Physics  and  Physics  Review  E.  More  details  can  be  found 
in  those  publications. 
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11  Galerkin  Method  for  Wave  Equations  with  Uncertain  Coefficients 


In  recent  years  there  is  a  growing  interests  in  studying  efficient  numerical  methods  for 
solving  differential  equations  with  random  inputs.  The  Polynomial  Chaos  (PC)  based 
methods  have  received  intensive  attention.  The  original  PC  method  was  developed  by 
R.  Ghanem,  cf.  [11],  and  was  inspired  by  the  Wiener  chaos  expansion  which  uses  Her- 
mite  polynomials  of  Gaussian  random  variables  to  represent  random  processes  [31]. 
Later  the  approach  was  extended  to  generalized  Polynomial  Chaos  (gPC)  where  gen¬ 
eral  orthogonal  polynomials  are  adopted  for  improved  representations  of  more  general 
random  processes  [32],  With  PC/gPC  serving  as  a  complete  basis  to  represent  random 
processes,  a  stochastic  Galerkin  projection  can  be  used  to  transform  the  (stochastic) 
governing  equations  to  a  set  of  deterministic  equations  that  can  be  readily  discretized 
via  standard  numerical  techniques.  Although  such  a  Galerkin  approach  is  effective 
in  many  problems,  cf.  [10,?],  its  application  to  hyperbolic  problems  has  been  limited 
as  of  now.  We  believe  that  the  primary  reason  is  that  the  properties  of  the  system 
of  equations  resulting  from  a  Galerkin  projection  is  not  fully  understood.  (When  the 
uncertainty  does  not  change  the  direction  of  the  characteristics,  the  Galerkin  system 
can  be  shown  to  be  hyperbolic  and  solved  in  a  straightforward  manner  [3].) 


We  discuss  in  this  paper  the  application  of  the  gPC  Galerkin  method  to  the  simula¬ 
tions  of  hyperbolic  systems  that  contain  uncertainties.  In  general  these  uncertainties 
may  enter  through  initial  conditions,  boundary  conditions  or  through  uncertainties 
in  the  coefficients  of  the  problem.  Here  we  deal  with  the  case  that  the  coefficients  are 
functions  of  random  variables.  In  particular  we  use  a  scalar  wave  equation  as  a  model 
and  study  the  situation  in  which  the  inflow-outflow  conditions  change  as  a  function 
of  a  random  variable.  The  problem  is  whether  it  is  possible  to  impose  boundary  con¬ 
ditions  on  the  deterministic  system,  consistent  with  the  boundary  conditions  of  the 
original  equation. 


We  show,  in  this  paper,  that  the  deterministic  system  is  a  symmetric  hyperbolic 
system  with  positive  as  well  as  negative  eigenvalues.  A  consistent  and  stable  method 
of  imposing  the  boundary  conditions  is  outlined.  The  boundary  conditions  are  not 
satisfied  exactly  at  the  boundaries  but  rather  to  the  order  of  the  scheme.  Convergence 
of  the  scheme  is  established. 


The  paper  is  organized  as  following.  In  Section  12  we  present  the  model  problem  of  a 
scalar  hyperbolic  equation  where  the  wave  speed  is  a  random  variable.  A  consistent 
set  of  boundary  conditions  are  presented  for  the  deterministic  system  resulted  from  a 
gPC  Galerkin  procedure,  and  we  prove  convergence  of  the  scheme.  In  Section  13  we 


62 


present  numerical  results  to  support  the  theory. 


12  Model  problem:  Scalar  Wave  Equation  with  Uncertainty 


A  simple  scalar  equation  that  illustrates  the  difficulties  in  applying  the  (generalized) 
Polynomial  Chaos  to  hyperbolic  equations  is: 


du(x,t,y ) 
dt 


c(y ) 


du(x,t,y ) 
dx 


x  G  (—1, 1),  t  >  0, 


(81) 


where  c(y)  is  a  random  transport  velocity  of  a  random  variable  y  e  12  in  a  properly 
defined  complete  random  space  with  event  space  12  and  probability  distribution  func¬ 
tion  p(y).  With  this  the  expectation  of  a  given  function  is  E [f(y)]  =  f  f(y)p(y)dy.  At 
this  stage  we  would  like  to  mention  that  we  can  consider  (81)  as  a  system  where  c  is 
a  symmetric  matrix  and  obtain  similar  results.  For  simplicity  we  stay  with  the  exam¬ 
ple  above  to  highlight  the  fundamental  properties.  The  physical  domain  is  bounded, 
(— 1, 1)  upon  proper  scaling,  so  that  we  can  study  the  effects  of  boundary  conditions. 


The  initial  condition  is  given  by 


u(x,0,y)  =  u0(x,y). 


(82) 


The  boundary  conditions  are  more  complicated  as  they  depend  on  the  sign  of  the 
random  transport  velocity  c(y).  A  well  posed  set  of  boundary  conditions  is  given  by: 


u{l,t,y)  =  uR{t,y),  c(y)>  0, 

u(-l,t,y)  =  uL(t,y),  c(y)  <  0. 


(83) 


Equations  (81)-(83)  complete  the  setup  of  the  problem. 
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12.1  Polynomial  Chaos  Galerkin  Approach 


Following  the  standard  gPC  expansion,  we  assume  that  u(x,  t,  y )  is  sufficiently  smooth 
in  y  and  has  a  converging  expansion  of  the  form 

OO 

u(x,  t,y)  =  t)Pk(y),  (84) 

k=0 

where  the  polynomials  Pk(y )  correspond  the  distribution  of  the  random  variable  y 
and  satisfy  the  following  orthogonality  relation 

E[PkPi]  =  j  Pk(y)Pi(y)p(y)dy  =  8kh  Vk,l,  (85) 

where  8ki  is  the  Kronecker  delta  function.  Note  the  polynomials  are  normalized.  The 
commonly  seen  correspondences  between  the  polynomials  Pk{y)  and  the  distribu¬ 
tion  of  the  random  variable  y  include  Hermite-Gaussian  (the  original  PC  expansion), 
Legendre-uniform,  Laguerre-Gamma,  etc.,  cf,  [32,?].  For  simplicity  we  will  discuss 
in  this  paper  the  case  of  random  variable  y  with  beta  distribution  in  (—1, 1)  (upon 
proper  scaling).  In  this  case  the  expansion  functions  Pk  are  the  (normalized)  Jacobi 
polynomials.  (Note  this  includes  the  special  case  of  Legendre  polynomials  with  uni¬ 
formly  distributed  random  variable  y.)  For  the  converged  series  (84),  we  also  assume 
that  the  expansion  coefficients  decay  fast  asymptotically,  i.e., 

ll%0M)lli  <  (86) 

where  K,  m  >  0  are  constants  and  the  ||  •  ||i  norm  is  defined  as 

||«jOM)||?  =  dx ■  (87) 

We  also  use  ||  •  || 2  to  denote  the  standard  L2  norm,  i.e.,  ||/(a;)|||  =  f\  f2{x)dx. 

By  utilizing  the  expansion  (84)  and  employing  a  Galerkin  projection,  it  is  straight¬ 
forward  to  verify  that  the  coefficients  v,j(x,t )  satisfy  the  following  infinite  system  of 
equations 
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d Uj(x,  t ) 


dt 


5Z  ax 


uk~ 


k= 0 
rl 


duk(x,t ) 
dx 


J  =  0,  •  •  • ,  oo 


aj,k  =  J  c(y)Pj(y)Pk(y)p(y)dy. 


(88) 

(89) 


The  equations  for  the  first  (N  +  1)  coefficients  can  be  written  as 


duj(x,t )  A  duk(x,t)  “  duk(x,t)  . 

=  2^aj,k — ^ - HZ — ,  J  =  0 


dt 


fc=0 


fc=W+l 


(9a; 


(90) 


In  the  gPC  Galerkin  method  we  seek  an  approximation  to  the  true  solution  via  a 
finite-term  gPC  expansion 


N 

v(x,t,y)  =  J2^k(x,t)Pk(y) 


k= 0 


(91) 


and  project 

dv(x,t,y)  dv(x,t,y) 

sr-  ~cM^^  =  0 

onto  the  subspace  spanned  by  the  first  (N  +  1)  gPC  basis  polynomials  and  obtain  the 
following  system 


dvj(x,  t)  _  A  dvk(x,t) 

~  ahk~ 


dt 


k= 0 


dx 


j  =  0, . . . ,  N, 


(92) 


where  ahk  are  defined  as  in  (89).  If  we  denote  by  A  the  ( N  +  1)  x  ( N  +  1)  matrix 
whose  entries  are  {ayfc}o<j,fc<iv  and  v  =  (h0,  •  •  •  ,  vn)T  a  vector  of  length  (N  + 1),  then 
system  (92)  can  be  written  as 


dv(x,t )  dv(x,t) 

dt  dx 


(93) 


Note  that  from  the  definition  aj  k  =  akj,  i.e.,  A  =  AT,  the  system  (93)  is  therefore 
symmetric  hyperbolic ,  this  is  consistent  with  the  fact  that  the  original  equation  (81) 
is  hyperbolic  for  each  realization  of  y. 
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12.2  Eigenvalues  of  the  PC  Galerkin  Equations 


A  less  trivial  question  is  the  nature  of  the  inflow-outflow  boundary  conditions.  The 
boundary  conditions  for  the  original  scalar  equation  (81)  depend  on  the  particular 
realization  of  the  random  variable  y  (see  (83)).  However  upon  the  Galerkin  projection 
in  the  random  dimension  the  deterministic  system  (93)  is  independent  of  y.  In  The¬ 
orem  1  we  investigate  how  the  inflow-outflow  conditions  are  reflected  in  the  system 
(93). 


Theorem  1: 


Consider  the  deterministic  system  (93)  where  the  coefficients  are  defined  in  (89).  Then 
if  c(y )  A  0  (reps.  c(y)  <  0)  for  all  y,  then  the  eigenvalues  of  A  are  all  non-negative 
(reps,  non-positive);  if  c(y)  changes  sign,  i.e.,  c(y)  >  0  for  some  y  and  c(y)  <  0  for 
some  other  y,  then  A  has  both  positive  and  negative  eigenvalues  for  sufficiently  large 
N. 


Proof: 


First  let  us  consider  the  case  of  c(y)  >  0. 


Let  (3(y )  be  a  random  variable  with  an  expansion  /3(y )  =  Jfk=0bkPk(y)-  Let  b  = 
(bo,  ■  ■  ■  ,  &tv)t  be  the  coefficient  vector  with  length  (N+l).  Note  here  b  is  an  arbitrary 
vector.  Then 


N  N 


bJ  Ab  =  ^2^2  bjajtkbk 

j= 0  fc= o 

=  J2  J2  bi  (/x  c(y)Pj(y)Pk(y)p(y)dy^j  bk 


j= 0  k= 0 

=  J_iP2(y)c(y)p(y)dy- 


Since  c(y)  is  non-negative 


brAb  >  0 


for  all  b,  thus  all  the  eigenvalues  of  A  are  non-negative.  The  case  of  c(y )  <  0,Vy 
follows  similarly. 

A  more  interesting  case  is  when  c(y)  changes  sign.  Let  us  devide  the  domain  hi  where 
y  belongs  into  the  following  non-overlapping  open  sets:  =  hi  U  hj  be  defined  as 

the  subdomain  of  y  where  c(y)  >  0  and  is  the  subdomain  of  y  in  which  c(y)  <  0. 
Let  us  also  define  7 (y)  be  a  smooth  function  such  that 


7  (y)>s,  yen  1,  (94) 

o<7  (y)<8,  yen2,  (95) 

i(y)  =  0,  yen3.  (96) 

Let  (3N(y)  be  the  best  polynomial  approximation  of  degree  N  to  \J^{y)  such  that 

max  \(3%(y)  —  7(2/) |  <  e,  (97) 


where  N  is  sufficiently  large  such  that 


c  fni  c(y)p(y)dy 
Jn\c(y)\p(y)dy 


(98) 


Then 


pN(y)c(y)p(y)dy=  7  {y)c(y)p(y)dy+  (P2N(y)-i(y))c(y)p(y)dy 

Jn  Jn+  Jn+ 

+  [  fd2N(y)c(y)p{y)dy 
Jn% 

>  7 (y)c(y)p(y)dy-  /  (P%{y)  -  j(y))c(y)p(y)dy 

Jn 1  Jn+ 

-  [  (d2N{y)c{y)p{y)dy  ■ 

Now 

/  fd2N{y)c{y)p(y)dy  <  e  /  \c(y)\p(y)dy 
J  n3  J  n3 

and 

/  (P%(y)  -  rr(y))c(y)p(y)dy  <  e  /  \c(y)\p(y)dy 
Jn+  Jn+ 

and  therefore 

[  fd2N{y)c{y)p{y)dy  >  o 
Jn 
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under  the  condition  (98).  Thus  there  exists  a  polynomial  ^n(v)  with  expansion  co¬ 
efficients  b  =  (bo,  ■  ■  ■  ,b]y)T  such  that  h1  Ab  is  positive.  Similarly,  if  c(y)  is  negative 
in  a  subinterval  there  exists  a  polynomial  / 3n(u )  with  sufficiently  large  N  such  that 
f\  /3jf(y)c(y)p(y)dy  is  negative.  Thus  the  matrix  A  has  positive  and  negative  eigen¬ 
values.  This  concludes  the  proof. 


12.3  Boundary  Conditions  and  Convergence 


We  now  turn  to  the  issue  of  imposing  the  boundary  conditions.  Since  A  is  symmetric 
there  is  an  orthogonal  matrix  S7  =  S_1  such  that 

StAS  =  A, 

where  A  is  a  diagonal  matrix  whose  entries  on  the  eigenvalues  of  A,  i.e., 

A  diag(A0, . . . ,  Aj+, . . . ,  A^_, . . . ,  A/y). 

Here  the  positive  eigenvalues  occupy  indices  j  =  0 , ,j+,  the  negatives  ones  j  = 
j_, . . . ,  N,  and  the  rest,  if  exist,  are  zeros.  Obviously,  j+,j-  <  N. 

Denote  by  q  =  (q0, ...,  qN)T  =  STv,  i.e., 


N 

Qj(x,t )  =  J2sk,jVk(x,t), 

k= 0 

where  Sj^  are  the  entries  for  S,  then  we  obtain 

aq^t)  =  Aaq M  (99) 

at  OX 

The  boundary  conditions  of  this  diagonal  system  are  determined  by  the  sign  of  the 
eigenvalues,  i.e.,  we  need  to  specify 


N 

&(M)  =  X)  skjUk(l,t),  J  =  0,  ...,j+, 
k= 0 

(100) 

N 

k= 0 
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Here  the  coefficients  Uk  at  the  boundaries  are  determined  by  the  exact  gPC  projection 
of  the  boundary  conditions  of  u,  i.e., 


UR{t,  y)  =  J2  “fc(1>  t)pk(y), 

3=0 

oo 

uL(t,  y)  =  t)Pk(y )• 

3=0 

Subsequently  the  boundary  conditions  for  the  gPC  Galerkin  system  of  equations  (93) 
are  specified  as 

v(l>  t)  —  Sq(l,  t)y  v(-l,t)  =  Sq(-l,t).  (101) 

Note  the  above  specification  of  boundary  conditions  via  (100)  and  (101)  implicitly 
satisfy  the  following  relation 


N 


N 


)  )  ■Sfcjhfc(l)  t)  ^  )  Sfcj'U/c(  1,  f), 

fc=0  k=0 


N  N 

)  )  s fcyhfc(  1;  t)  ^  )  ■Sfcyh/C(  1,  f), 


fc=0 


fc=0 


J  =  0, 

J  =]-,■■■,  N. 


For  vanishing  eigenvalues,  if  they  exist,  no  boundary  conditions  are  required. 


Theorem  2 


Consider  the  hyperbolic  equation  equation  (81)  where  y  is  a  random  variable  with 
beta  distribution  in  (—1,1).  Let  u(x,t,y )  be  the  solution  of  (81)  whose  exact  gPC 
expansion  is  (84)  and  let  v(x,t,y )  be  the  (N  +  l)-term  gPC  solution  (91)  solved  via 
the  Galerkin  system  (93)  with  boundary  conditions  given  in  (101).  Then  for  any  finite 
time  t 


E 


u  —  v 


2 

2 


Vj(x,t))2dx'Sj  < 


K 


yy2m— 1 


t. 


(102) 


Note  the  linear  growth  in  time. 


Proof 
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Let 


ej(x,  t )  =  Uj(x,  t)  —  Vj(x,  t ),  j  —  0, . . . ,  N. 
^Frorn  (90)  and  (92)  we  have 


dej(x,t)  A  dek(x,t) 

=  1 a3,k - - - h 


dt 


k= 0 


dx 


aj,k 

k=N+ 1 


dx 


J  =  0, 


,  N. 


(103) 


Denote  by  e  =  (eo,  •••,  ejv)T  and  let  d  =  STe,  then  we  obtain 

dd  4  dd  „ 

~dt~Ad^+  * 

where  the  residual  vector  R  =  (R0,  ■  ■  ■  ,  Rn)t  is 


N  oo 

Rj(x,t)  =  j2  suaw 

1=0  k=N+ 1 


dx 


(104) 


(105) 


By  multiplying  (104)  by  dT  and  integrating  in  x  one  gets 


1  d 

2  dt 


dTd  dx 


/- 1 


i  N  ! 

9  Ai  (C?j(1>  *)  -  dj(- !,  *))  +  /  dTRoh 
j=o 


(106) 


^ From  the  boundary  conditions  (101)  it  follows  that  if  Aj  >  0  then  dj(l,t)  =  0  and 
if  A  j  <  0  then  dj(—l,t)  =  0.  Thus  the  first  term  in  the  right  hand  side  of  the  above 
equation  is  negative.  This  leads  to 

\li !-idTddx=  ISl|d|'2  “  l|d|1 ' l|R|1, 


where 


and  ||d||  is  defined  similarly. 


Rj(x,  t)dx 


Thus 


dt 


|d[|  <  IIRI 


70 


and 


||d(a;,  t)  II  <  max  ||R||  •  t 

Since  the  matrix  S  is  unitary  and  the  elements  of  the  matrix  A  are  bounded  then 

OO 

l|Rf<  E  INI?  (107) 

j=N+ 1 

and  the  proof  is  established  under  assumption  (86). 


13  Numerical  Results 


In  this  section  we  present  a  few  numerical  examples  to  support  the  theoretical  results 
derived  above.  In  all  of  the  following  computations,  we  have  used  sufficiently  foie 
resolutions  in  physical  space  and  time  domain,  such  that  the  spatial  and  temporal 
errors  are  negligible.  In  all  computations,  y  is  a  random  variable  uniformly  distributed 
in  (—1,1)  and  thus  Pk(y)  are  (normalized)  Legendre  polynomials. 


13.1  Periodic  Problem 


We  first  consider  problem  (81)  with  a  periodic  boundary  condition  in  physical  space. 
Subsequently  the  gPC  Galerkin  system  (93)  requires  periodic  boundary  conditions 
that  can  be  trivially  implemented.  Therefore  no  errors  will  be  induced  by  specifying 
boundary  conditions  via  (101).  Let  us  consider 


ut(x,  t,  y)  =  yux(x,t,y),  0  <  x  <  2n,  t  >  0, 
u(x,  0,y)  =  cos(x),  0  <  x  <  2i r, 


(108) 


The  exact  solution  is  uex  =  cos (x  —  yt).  In  Figure  25  we  plot  the  evolution  of  the  mean 
square  solution  E[||tt|||]  =  /  f\  u2(x,  t,  y)p{y)dydx  and  its  numerical  solution  via  gPC 
Galerkin  method.  We  observe  that  there  is  a  finite  time  where  the  numerical  solutions 
lose  accuracy,  i.e.,  errors  become  0(1).  The  size  of  the  time  domain  in  which  the  errors 
remain  small  grows  almost  linearly  as  the  orders  of  gPC  expansion  are  increased. 
This  observation  can  be  cross-examined  by  comparing  the  mean-square  errors  at 
different  time  level,  as  shown  in  Figure  26.  We  observe  that  with  sufficiently  high 
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orders  of  gPC  expansions,  exponential  error  convergence  can  be  achieved.  However, 
as  time  increases,  the  critical  orders  of  expansions,  beyond  which  errors  start  to 
decay  exponentially  fast,  increase  linearly.  All  of  these  results  support  the  convergence 
analysis  (102)  where  a  linear  error  growth  in  time  exists. 


Fig.  25.  Evolution  of  errors  in  mean-square  norm  over  time. 


13.2  Boundary  Conditions  and  Discontinuity  in  Random  Space 


We  now  study  a  wave  equation  with  a  random  wave  speed  that  changes  signs  and 
also  contains  a  discontinuity  in  the  random  space 


c(y)'U3;, 
u(x,  0,  y)  =  sin  (/ex), 
u(x,0,y)  =  sin(2/cx), 


— 1  <  x  <  1,  t  >  0, 
— 1  <  x  <  1,  y  >  0, 
— 1  <  x  <  1,  y  <  0. 


(109) 
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Fig.  26.  Convergence  of  mean-square  errors  with  increasing  orders. 

Here  c(y)  =  cry  with  0  <  a  <  1  controlling  the  variability  of  the  random  input  and 
k  >  0  is  a  real  constant.  We  prescribe  boundary  conditions  as 

u{l,t,y)  =  sin  |/c  (1  +  c(y)t)], 
u(- 1,  t,  y)  =  sin[2/c(— 1  +  c(y)t)\, 

The  exact  solution  of  (109)-(110)  is  ue(x,t,y )  =  sin[fi;(a;  +  c(y)t )]  for  y  >  0  and 
ue(x,t,y)  =  sin [2 n(x  +  c(y)t)\  for  y  <  0.  Note  the  solution  is  discontinuous  in  term  of 
?/,  although  each  realization  of  y  is  a  smooth  function  in  x. 


V  >  0, 

y  <  o. 


(110) 


The  numerical  solutions  are  solved  with  a  =  0.5  and  k—1.  The  numerical  boundary 
conditions  are  implemented  via  the  eigenvalue  analysis  (101).  For  numerical  solu¬ 
tions  of  gPC  order  N,  we  examine  three  error  measures:  error  in  mean  emeSLn(N,t)  = 
maxx  |E(u)  —  E(we)|,  error  in  standard  deviation  (STD)  estd  (N,t)  =  m&xx  \  av  — 
<jUe  | ,  and  the  mean-square  error  (IV,  i)  =  maxa,(E[(n  —  we)2])1//2.  Numerical  sim¬ 
ulations  are  conducted  up  to  t  =  1,  and  we  define  convergence  rate  as  r(N )  = 
[ln(e(iV))  —  ln(e(M))]  /  [In (AT)  —  ln(M)]  for  expansion  orders  N  >  M  >  1,  where  e  is 
one  of  the  three  error  measures. 


Figure  27  shows  the  convergence  of  the  three  errors  with  increasing  order  of  Legendre 
expansions.  In  this  case,  we  observe  different  convergence  properties  between  even  and 
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odd  orders  of  expansions,  although  they  appear  to  have  similar  asymptotic  conver¬ 
gence  rate.  Note  that  such  different  error  behaviors  between  even  and  odd  expansions 
can  be  seen  in  classical  spectral  methods,  cf.  [14].  The  errors,  along  with  their  con¬ 
vergence  rates,  are  tabulated  in  Table  XI  and  Table  XII,  for  odd  and  even  orders  of 
expansions,  respectively.  No  exponential  convergence  is  achieved,  as  opposed  to  that 
in  the  earlier  examples.  Also,  the  weak  error  measures  (error  in  mean  and  error  in 
STD)  converge  more  rapidly  than  the  strong  error  measure  in  term  of  mean-square. 


Fig.  27.  Errors  for  odd  and  even  orders  of  Legendre-chaos  expansions. 


The  slower  convergence  is  due  to  the  discontinuity  in  random  space,  and  is  manifested 
in  Fig.  28,  where  the  numerical  solution  v(x,t,y)  is  shown  at  location  x  =  0.454  with 
N  =  21  order  of  expansion.  Fig.  28(a)  shows  the  approximation  at  t  —  0,  i.e.,  the 
initial  condition,  and  Fig.  28(b)  shows  the  numerical  solution  at  t  —  1.  The  Gibbs’ 
oscillations  around  the  discontinuity  at  y  —  0  are  clearly  visible. 


14  Summary 


The  properties  of  (generalized)  Polynomial  Chaos  method  for  uncertainty  analysis  of 
hyperbolic  equations  are  studied.  We  show,  via  a  simple  model  problem  of  a  scalar 
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Table  XI 

Errors  for  odd-order  Legendre-chaos  expansions  and  their  convergence  rate. 


Table  XII 

Errors  for  even-order  Legendre-chaos  expansions  and  their  convergence  rate. 


wave  equation  with  random  wave  speed,  some  prominent  features  of  the  resulting  de¬ 
terministic  system  of  equations  obtained  by  a  Galerkin  projection  in  random  space. 
We  proved  the  existence  of  both  positive  and  negative  eigenvalues  when  the  wave 
speed  changes  sign  in  random  space  and  presented  a  consistent  and  stable  method 
for  imposing  boundary  conditions  for  the  deterministic  equations.  The  gPC  Galerkin 
method,  with  the  proper  boundary  treatment,  is  shown  to  be  convergent.  Further¬ 
more,  the  error  contains  a  linear  growth  in  time  which  is  independent  of  the  bound¬ 
ary  conditions.  We  remark  that  although  the  linear  wave  equation  considered  here 
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Fig.  28.  Numerical  approximations  of  v,N(x,t,y)  at  x  =  0.454  with  N  =  21.  (a)  t  =  0;  (b) 
t  =  1. 

is  rather  simple,  it  possesses  one  of  the  key  issues  in  applying  gPC  Galerkin  method 
to  hyperbolic  problems  -  the  proper  way  to  enforce  boundary  conditions  when  the 
characteristic  wave  changes  directions  in  random  space.  This  issue  is  addressed  here 
and  it  opens  up  the  possibility  of  applying  gPC  Galerkin  method  to  other  hyperbolic 
problems  with  uncertainty,  e.g.,  nonlinear  wave  equations,  Maxwell  equations,  etc. 
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